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ABSTRACT 


The objective of this study was to develop a high-resolution-explicit-multi-block 
numerical algorithm, suitable for efficient computation of the three-dimensional, time- 
dependent Euler and Navier-Stokes equations. The resulting algorithm has employed a 
finite volume approach, using MUSCL-type differencing to obtain state variables at cell 
interface. Variable interpolations were written in the K-scheme formulation. Invtscid 
fluxes were calculated via Roe’s flux-difference splitting, and van Leer’s flux-vector 
splitting techniques, which are considered state of the art. The viscous terms were 
discretized using a second-order, central-difference operator. 

Two classes of explicit time integration has been investigated for solving the com- 
pressible inviscid/viscous flow problems —two-stage predictor-corrector schemes, and 
multistage time-stepping schemes. The coefficients of the multistage time-stepping 
schemes have been modified successfully to achieve better performance with upwind 
differencing. A technique was developed to optimize the coefficients for good high- 
frequency damping at relatively high CFL numbers. Local time-stepping, implicit resid- 
ual smoothing, and multigrid procedure were added to the explicit time stepping scheme 
to accelerate convergence to steady-state. The developed algorithm was implemented 
successfully in a multi-block code, which provides complete topological and geometric 
flexibility. The only requirement is C° continuity of the grid across the block interface. 

The algorithm has been validated on a diverse set of three-dimensional test cases of 
increasing complexity. The cases studied were: (1) supersonic comer flow; (2) supersonic 
plume flow; (3) laminar and turbulent flow over a flat plate; (4) transonic flow over an 


ONERA M6 wing, and (5) unsteady flow of a compressible jet impinging on a ground 
plane (with and without cross flow). The emphasis of the test cases was validation of 
code, and assessment of performance, as well as demonstration of flexibility. 
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CHAPTER 1 
INTRODUCTION 


1.1 Historical Background 

Recent developments in numerical methods and their applications permit the solving 
of complex, realistic geometries and configurations for compressible flows. Currently, 
computational fluid dynamics (CFD) is used effectively to unravel and elucidate fluid 
flow phenomena which are difficult to study in the laboratory. The demand to solve 
finely detailed models of physics has challenged many researchers to come up with new 
and efficient tools. This demand has resulted in revolutionary concepts in computer 
architecture designs and software development. 

The birth of CFD can perhaps be linked to the early work of the English mathe- 
matician Richardson in 1917 [1]. He attempted to integrate the meteorological equations 
numerically. It is interesting to note that he started this process, which evolved into a 
new science, as an ambulance driver during World War I. He made the computations 
by hand, [1]. His attempts were unsuccessful due to a limited theoretical understanding 
of the stability of numerical methods, and to a lack of computing power. Richardson’s 
failure; outlined the areas which needed to be developed. In 1928 Courant, Friedrichs, 
and Lewy [2] introduced their famous stability condition, which became subsequendy 
the CFL number, and represented a landmark mathematical result that has had a massive 
impact on computational research. 

The practical birth of CFD came in the late 1960’s when significant computing power 
became available. Since then, there has been considerable progress in the field of CFD. 


1 


The growing field of aerodynamics, and the aviation industry have been the catalyst for 
the revolutionary force of CFD. In this section, a brief review of previous Computational 
Fluid Dynamics work related to the present work is presented. For a more general review 
of CFD, the interested reader should review Refs. [1, 3-6]. 

One of the first major advances in Computational Fluid Dynamics was the work by 
Hess and Smith [7], who introduced panel methods to solve the linearized potential flow 
equation. Later the panel method was extended to lifting flows by Rubbert and Saaris 
[8] and supersonic flows by Woodward [9]. In 1986, Kandil and Yates [10] extended 
the method to solve the steady, full-potential equation. In 1987, Kandil and Hong [11], 
successfully formulated the vortex-panel method in a moving frame of reference. 

In the early seventies, two major breakthroughs were reported which allowed the 
solution of non-linear mathematical models. Murman and Cole [12], devised the idea of 
mixed differencing (central differencing in subsonic regions, and forward or backward 
differencing in supersonic regions of flow). They employed a line relaxation method for 
the entire flow field, which was partly elliptic and partly hyperbolic. Their work, and 
the work of Jameson [13], was the catalyst for developing two- and three-dimensional 
algorithms using the Small Disturbance Equation, and the Full Potential Equation. 
An interesting review of the memoirs of Murman and Cole is presented in a review 
paper by Hall [14], The second major breakthrough was the work by Magnus and 
Yoshihara [15]. They advanced the Euler Equation in time towards a steady-state, 
thus transforming a mixed elliptic- hyperbolic problem into a purely hyperbolic one. 
Another landmark in the history of CFD came in 1970, when McCormack introduced 
his widely used predictor-corrector explicit difference scheme [16]. Subsequently, in 
1981, McCormack [17] developed an implicit analogue of his explicit finite difference 
method. In 1975, Wanning and Beam [18] introduced a fully upwind predictor-conector 


2 



method, which is similar to the McCormack method. Briley and McDonald [19], and 
Beam and Warming [20, 21] employed an Alternating Direction Implicit (ADI) scheme 
for solving the Euler and Navier-Stokes equations. The roots of ADI schemes trace back 
to Peaceman and Rachford [22], Douglas [23], and Douglas and Gunn [24], Steger [25], 
adapted the Beam and Warming scheme to general curvilinear coordinates. ADI evolved 
to an effective tool and currently is employed in state-of-the-art codes designated ARC2D 
and ARC3D [26]. 

On the other hand, another important family of time integration schemes —explicit, 
multistage time-stepping schemes (Runge-Kutta methods) — started to evolve in the early 
eighties. Jameson, Schmidt, and Turkel [27], introduced explicit, multistage time-stepping 
schemes, to the CFD community. Explicit-multistage schemes were developed further, 
and have been applied successfully to compute solutions to the Euler, and Navier-Stokes 
equations, for two- and three-dimensional problems [28-35]. Explicit schemes combine 
naturally with accelerating techniques such as: local time-stepping, residual smoothing, 
and multigrid accelerating techniques. They are also well suited for parallel computing 
[36, 37]. 

The restriction on the time step for explicit schemes was the catalyst to develop 
implicit schemes. Implicit schemes require more computation per time step (iteration), 
but allow a larger time step to be used. The implicit time integration scheme may be 
stable for any step size, according to linear theory, yet it is limited in practice by the non- 
linearity of the governing equations. Due to simplifications made during the development 
of these methods (linearization) and the frequent use of explicit boundary conditions, the 
maximum allowable Courant-Friedrichs-Lewy number ( CFL ) is reduced. To date, the 
relative merits of implicit and explicit schemes are still an open debate for steady and 
unsteady flow calculations. 
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Implicit residual smoothing extends the stability limit, and improves the damping 
properties of the multistage time-stepping scheme. Lerat [38], introduced the idea 
of residual smoothing for the Lax-Wendroff scheme [39], Jameson and Baker [29] 
applied the idea of implicit residual smoothing in conjunction with the modified Runge- 
Kutta schemes. This procedure was developed further in Refs. [28, 33, 40-42], where 
they employed a central-implicit-residual-smoothing operator. The use of an upwind- 
residual-smoothing operator was employed by von Lavante and Gronner [43], and Blazek 
et al. [44], 

Multigrid acceleration techniques were developed originally by Fedorenko [45, 46] 
starting in 1961. Subsequently Brandt [47] applied the technique to an elliptic set of 
equations. The work by Brandt and many others has led to the popular use of multigrid 
by many in the fields of applied mathematics and computational engineering. Excellent 
developments of the multigrid technique can be found in Refs. [48-50]. Multigrid 
was used successfully for solving the potential, Euler, and Navier-Stokes equations. 
Refs. [51-55]. Multigrid acceleration techniques performed well when combined with 
central-difference methods, but the convergence rate deteriorated with upwind spatial 
operators because they are less dissipative. One must ensure that the basic upwind 
scheme exhibits good damping of high frequencies on both fine and coarse meshes. An 
attempt to derive a mathematical operator to eliminate the high-frequency components 
of the error should be pursued. 

In the early 1980’s, computers were powerful enough to permit the computation of 
solutions to the Euler equations. A new wave of inviscid upwind and central-difference 
schemes evolved. Upwind schemes attempted to construct the flux by modelling the un- 
derlying physics, as dictated by the sign of characteristic waves, while central-difference 
schemes computed the interface flux as an average of the two adjacent cells, disregarding 
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characteristic theory. A comparison between central-difference schemes, and upwind 
schemes and how they are related is given by Swanson and Turkel [56], 

A prevalent way to introduce upwinding into the governing systems for hyperbolic 
conservation equations has been to split the flux according to the characteristic speed 
(q, q±c)). Steger and Warming [57], were the first to devise a conservative-second- 
order-flux-vector splitting-upwind scheme, without the use of limiters, for the solution 
of the governing equations of gas dynamics. Anderson, Thomas, and van Leer [58], 
developed the Monotone Upstream-centered Scheme for Conservative Laws (MUSCL) 
approach with limiters which was incorporated in the Steger- Warming scheme. The 
MUSCL approach resulted in a better shock capturing capability. The main disadvantage 
of the Steger- Warming-flux-vector splitting scheme was that the backward and forward 
fluxes were not differentiable. This leads to oscillation at shocks, van Leer [59], devised 
an alternative splitting scheme. The advantage of van Leer’s flux-vector splitting over 
the Steger- Warming flux-vector splitting scheme, was that the split flux-vectors were 
smooth and had smooth first derivatives with respect to the Mach number, so that their 
eigenvalues were also smooth [58]. 

The inviscid flux can be split in a number of ways. The Split Coefficient Matrix, 
(SCM) as introduced by Chakravarthy, Anderson, and Salas [60], is a natural way of 
splitting the flux based on the sign of the eigenvalues of the governing system of 
equations. A similar scheme that is based on the theory of characteristic is Morreti’s 
A-scheme [61]. Both the SCM- and A-schemes have been applied to the non-conservative 
form of the governing equations, and require shock-fitting techniques in the presence of 
shocks. The conservative form of the governing equations permit shock waves to be 
captured as weak solutions to the governing equation [39, 62, 63], thus avoiding the 
difficulty of applying shock-fitting techniques. 
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In 1959, Godunov [64], introduced the idea of advancing in time by solving the 
Riemann problem at each cell. This technique has been extended to higher order schemes 
which are known today as Gudonov-type schemes, [65-71]. A review of Gudonov-type 
schemes is presented by Roe [5], and Yee [72], Currently, upwind schemes are being used 
on a regular basis for computing solutions to the Euler and Navier-Stokes equations. They 
have been implemented and validated in several state-of-the-art codes, such as CFL3D 
[73], ISAAC [74], PAB3D [75 ] and FTNS3D [76] . 

Alternatively, Jameson, Schmidt, and Turkel [27] have introduced multistage time- 
stepping schemes, coupled with a central-difference operator and explicitly added dissi- 
pation terms. The explicit dissipation term was a blend of second-order-difference and 
fourth-order-difference terms. Second-order-difference terms suppress oscillations in the 
neighborhood of shock waves, while fourth-order-difference terms are crucial for the sta- 
bility and convergence to steady-state. Dissipation terms have been scaled by user defined 
coefficients. Detailed discussion of the influence of the dissipation terms on the perfor- 
mance and quality of steady-state solutions can be found in Kandil and Chuang [77], 
Rizzi [78], Pulliam [79], and Swanson and Turkel [51]. 

Currently, the state-of-the-art in computational fluid dynamics replaces scalar dissipa- 
tion with a matrix-valued dissipation function. Employing matrix dissipation enhances the 
shock capturing capabilities of the central-difference technique, and reduces the smearing 
of shocks and contact discontinuities which were characteristic of the original central- 
difference schemes [51]. Central -difference operators, coupled with a matrix-valued dis- 
sipation function, are nearly as accurate as upwind schemes, and have the merit of being 
computationally cheaper and easier to program [56]. 

The numerical dissipation terms play an important role in the success of the compu- 
tations by central-difference methods. For every new configuration, the exact (optimum) 
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level of artificial dissipation is not known a priori. The level of numerical dissipation 
can be turned up, by a novice user, to the point of masking the physics of the prob- 
lem. A certain level of expertise with central-difference schemes and with the physical 
problem of interest is required to select the optimum (acceptable) level of dissipation. 
Central-difference schemes have been applied in state-of-the-art schemes, TLNS3D [28], 
ARC2D and ARC3D [26] and FLOMG [32], 

The application of the above numerical methods to realistic three-dimensional con- 
figurations of significant geometric complexity is virtually impossible without the use 
of Domain Decomposition techniques. Here, the computational domain is divided into 
multiple blocks (zones) and the grid for each block is then generated. A computational 
grid of this type adapts more easily to the shape of the body as well as to the flow 
features. Typically, the transfer of information between the blocks is carried out explic- 
itly by ensuring the conservation of fluxes across the block interfaces. The consequence 
of this procedure, for an implicit operator, is a significant reduction in the maximum 
allowable CFL number. 

Generating a single body fitted grid for complex, three-dimensional realistic ge- 
ometries is a difficult task to perform; for some configurations it is almost impossible 
[80-82]. Several grid methodologies such as overlaid grids [83], patched grids [84], 
blocked grids [85], and unstructured grids can be applied to simplify the grid generation, 
provide geometric flexibility, and even provide mesh refinement. Several methods have 
been investigated for unstructured grids. Refs. [86-89]. These methods require more 
memory and computational time and fall short of their structured counterpart in terms of 
efficiency and accuracy [87]. The theory and algorithms for unstructured grids have to 
evolve before they can be used for solving practical three-dimensional problems. 
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Explicit-upwind schemes appear to be a good compromise between explicit-central- 
difference schemes, and implicit-upwind schemes. Schemes constructed along these lines 
combine the advantages of; simplicity; prudent use of computational resources, and 
accuracy in resolving the flow field. Upwind schemes are more complex and are usually 
reported to be better suited for compressible viscous computation. Upwind schemes are 
very effective in converging to steady-state on single grids of modest complexity. Most 
of the currently used upwind schemes are implicit. Explicit schemes require less memory, 
and are easily implemented in a multi-block environment. They are also naturally suited 
for implementation on massively parallel computer architecturs. The main drawback of 
explicit schemes is the limitation on the allowable time step. 

If the explicit time-stepping scheme is augmented with suitable accelerating tech- 
niques, such as local time-stepping, residual smoothing and multigrid acceleration, the 
explicit method will be superior to its implicit counterpart. Variable coefficient resid- 
ual smoothing will increase the stability range of the scheme, thus allowing the use of 
a higher CFL number (larger time step), which enhances the rate of convergence and 
removes the diffusion limit on the time step. Multigrid acceleration techniques will ac- 
celerate the convergence to steady-state by using large time steps on coarser grids, and 
help achieve convergence rates that are independent of the number of grid points [47] 

Recent advances in computer architecture and algorithmic tools open the door for 
a new wave of opportunities for constructing explicit, upwind-higher-order schemes. 
Currently the existence of robust, multi-block, explicit, upwind schemes that can be 
applied on a routine basis are not available. Upwind-high-order schemes are essential 
tools, required to capture complicated physical phenomena associated with problems of 
current interest. 
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Explicit-upwind-schemes are still in their infancy and many basic issues are yet to be 
settled. In order to lay the foundation for future research, a joint analytical and numerical 
study should be conducted to validate and demonstrate their capabilities and performance. 

1.2 Objective of Present Work 

The goal of the present work was to develop a general state-of-the-art, multi-block 
algorithm, capable of solving the governing equations of fluid motion efficiently, for 
a wide range of configurations with both internal and external flow. The requisite 
algorithm should be simple, efficient, and robust. It is required to damp the high frequency 
component of the error (necessary for multigrid) effectively, while acquiring low levels 
of numerical dissipation for accurate predictions of viscous effects, and still maintaining 
high resolution on stretched grids. The developed algorithm will be used subsequently 
to simulate complex three-dimensional, steady and unsteady flow problems. 

Hence, a control- volume, explicit-multistage-high-resolution upwind scheme, suitable 
for efficient computations using block structured grids, was desired. Upwind schemes 
were selected due to their high degree of reliability in viscous flow computations and 
their superior shock capturing capabilities [90]. The state variables at the cell interface 
have been determined by MUSCL interpolation using the so-called n scheme. Two state- 
of-the-art, upwind schemes: Roe’s flux-differencing and van Leer’s flux-vector splitting 
schemes, were u tiliz ed to evaluate the inviscid flux at the cell interface. The viscous 
stress and heat flux terms in the governing equations have been centrally differenced. 

In this study, the objective was to devise explicit, upwind time-stepping schemes 
that can be combined successfully with upwind-spatial operators. Explicit schemes 
have the merit of being computationally cheaper and easier to program and implement 
in a multi-block code. Two classes of upwind schemes: multistage time-stepping 
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schemes, and predictor-corrector schemes, were suggested and have been implemented 
in the developed algorithm. Modified Runge-Kutta methods with standard coefficients 
have been successful with central-difference spatial discretization. Yet, they have not 
performed as well with upwind differencing. The standard coefficients have to be 
modified to achieve better performance with upwind differencing. 

The next objective was to augment the explicit time-stepping schemes with accel- 
erating techniques, such as local time-stepping, implicit residual smoothing and the full 
approximation storage (FAS), to enhance the rate of convergence to steady-state. 

Current aerodynamics designs are often quite complex (geometrically). Flexible 
computational tools are needed for the analysis of a wide range of configurations with 
both internal and external flows. Hence, another objective was the implementation of the 
developed algorithm in a multi-block code to allow for greater geometric flexibility. 

The final goal was to validate the developed computer code on several test cases of 
interest to demonstrate and assess the predictive capability of the algorithm. The test 
cases considered were: corner flow, plume flow, laminar and turbulent flow over a flat 
plate, an ONERA M6 wing, and the unsteady three-dimensional flow of a jet impinging 
on a ground plane. 

1.3 Thesis Outline 

In chapter two, the mathematical formulation of the governing set of equations of 
motion (Reynolds -Averaged Navier-Stokes equation, Navier-Stokes equation, and Euler 
equation) are presented and discussed. Details of implementing the Baldwin-Lomax 
algebraic eddy viscosity turbulence model in the algorithm are presented. In chapter three, 
the finite volume formulation of the governing equations is presented. The MUSCL type 
differencing, and the type of discretization for the inviscid and viscous flux is discussed. 
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Two different upwind flux formulations —Roe’s flux-difference splitting, and van Leer’s 
flux-vector splitting — are presented, and practical issues concerning their implementation 
are discussed. In chapter four, the temporal discretization of the governing equations, 
which represents a major part of this work is presented. Two classes of explicit 
time integration schemes —multistage time-stepping schemes and predictor-corrector 
schemes— are presented and discussed. Details of optimizing the multistage explicit 
time integration scheme through local Fourier analysis of the scalar advection equation 
are presented. Accelerating techniques, including local time-stepping, residual smoothing, 
and multigrid acceleration techniques are presented in chapter five. In chapter six, the 
multi-block capability of the developed algorithm is presented. The interaction between 
multigrid and multi-block implementations are discussed. The boundary conditions 
employed, in the developed algorithm, are presented within the framework of multi-block. 
Several test cases of general interest to the computational fluid dynamics community were 
conducted to validate, demonstrate and assess the performance and predictive capability 
of the present algorithm. Results of these computations are reported in chapter seven. 
The Conclusions for the present research work, and recommendations for future research 
are presented in chapter eight. 
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CHAPTER 2 

GOVERNING EQUATIONS 


The governing equations were derived from the basic principles of conservation of 
mass, conservation of momentum, and energy. The conservation laws were then coupled 
with the thermodynamic properties and constitutive equations to yield the governing set 
of equations for fluid motion. The derivation of the governing equations can be found 
in [91, 92]. Three different sets of governing equations have been used pertaining to 
the different test cases investigated in this study. These sets of equations are the Euler 
equations, Navier-Stokes equations, and Reynolds-Averaged Navier-Stokes equations. 
Each set is represented ultimately as an algebraic set of equations. The three sets of 
governing equations have been implemented in the numerical algorithm. Coupled with the 
appropriate set of boundary conditions, the developed algorithm is capable of computing 
inviscid, laminar, and turbulent fluid flows numerically. 

2.1 Navier-Stokes Equations 

The time-dependent, compressible, three-dimensional, Navier-Stokes equations in 
Cartesian Coordinates, written in strong conservation form (neglecting the body forces 
and external heat sources) are: 
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where Q is the state vector of the dependent variables, given by: 


p 


' Density 

p u 


x — momentum 

pv 

► = i 

y — momentum 

pw 


z — momentum 

k E t 


Total energy , 


F, F v , G , G v , H, and H v are the dimensional flux-vectors. They are function of the state 
variable vector, Q, and are given by: 


F = { 


G= { 


H={ 


pu 1 + p 

puv 

pu w 

[e + p^U J 

V 

puv 
pv 2 + p 

pvw 

(E + pjd J 

w 

pu w 
pvw 

pw 2 + p I 

(s + pjw J 


F„ = { 


G v — \ 


X z ^ 


T~ — 

* y 

’ps ✓‘s. 

jy 

T'"'" 

V x 

+ V T — + iG — q 

1 tf M 1 II T * 


U T- 


x y 


y y 


y x 


y 7 


H v = < 


0 


X z 

'p*. 

JF z 


l UT~~+VT~~+WT~~-q: 


(2.3) 


(2.4) 


(2.5) 


The first row in the vector differential equation, eq. 2.1, is the conserved form of the 
continuity equation, while the fifth row is the conserved form of the energy equation. The 
second, third and fourth row are the conserved forms of the momentum equation in the x, 
y, and z directions, respectively. It should be emphasized that while the conservation of 
mass and energy are scalar equations, the conservation of momentum is a vector equation 
with three components. In the absolute sense only the x-momentum, y-momentum, and 
z-momentum are the Navier-Stokes equations of fluid mechanics. It is customary within 
the computational fluid mechanics arena to refer to eq. 2.1 as the Navier-Stokes equations. 
This terminology will be adopted in this study. 
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Fluid density is designated as p\ u, v, and w are the velocity components in the x, y, 
and z directions, respectively; E is the total energy; r f> are the components of the shear 
stress tensor; q x , q y , and q z are the components of the heat flux-vector in the x, y, and z 
directions respectively; and T is the temperature. The superscript in vector eq. 2.1 
indicates a dimensional quantity. 


In this study we assume that the stress is linearly dependent on the rate of strain; 
i.e., the Newtonian fluid assumption is adopted. The components of the viscous stress 
tensor in Cartesian coordinates are given by 
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( 2 . 6 ) 


y, is the first coefficient of viscosity, and A is the second coefficient of viscosity. To 
date, the value for A for air, and how to model it, especially for compressible flows, is 
still disputed [93]. In this study we employ Stokes hypothesis; i.e., we assume the bulk 


viscosity, K, is zero or negligible 


K = A + ^/i = 0 (2-7) 

O 

Stokes hypothesis is not necessarily endorsed, but for lack of a better model, it has been 
employed. It is understood that this assumption is not valid in shock regions and in regions 
of high gradients [94]. Invoking Stokes hypothesis yields the following expression for 
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the shear stress terms 
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The first coefficient of dynamic viscosity varies with temperature and can be approximated 
by an approximation of Sutherland’s law [93, 95] 



where fi ref and f re/ are the dynamic viscosity and temperature at reference conditions. 
This formulation is simple and gives reasonably good results. The parameter n is taken 


to be 0.76. 

The heat flux is modeled by Fourier’s law of heat conduction, where 
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where k is the thermal conductivity which will be represented by 
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Pr 

Here, Pr is the Prandtl number, and c p is the specific heat at constant pressure. The 
Prandd number is nearly constant for most gases (Pr = 0.72 for air). 

The equation of state for a perfect gas relates the pressure to density and temperature 
and is given by p = pRT, where H is the gas constant, which relates the specific 


heats for an ideal gas by: 



( 2 . 12 ) 
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with Cv representing the specific heat at constant volume. If we assume further a 
calorically perfect gas 



(2.13) 


where e is the internal energy. Now, neglecting the potential energy, the total energy, E, 
can be defined as the sum of the internal energy and kinetic energy; 


E = p 
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If we define 7 = as the ratio of specific heats, the above equations can be combined 

Cv 

to yield the relation between pressure and total energy 
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2.2 Normalization of the Governing Equations. 


Computing in an appropriate non-dimensional or normalized domain has the advan- 
tage that all variables are of the same order of magnitude, which enhances the performance 
and accuracy of numerical algorithms. Normalization eliminates the physical dimensions 
from the equations. Thus allowing general characteristic parameters such as Reynolds 
number, Mach number and Prandtl number to be changed independently. Hence para- 
metric studies on any of these characteristic parameters can be performed easily. 
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Different variables or combinations of variables can be used in the normalization 
procedure. In this study we define the non-dimensional flow variables to be 
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Here a re f is the reference speed of sound; H is the total enthalpy, and L is a reference 
length. The subscript ref indicates the reference condition. By substituting the non- 
dimensional flow variables into the Navier-Stokes equations, eq. 2.1 we get; 
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where Q, F, and F v are given by 
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where M re f, and R re f are the reference Mach number and Reynolds number respectively, 
while q Te f is the velocity magnitude at reference conditions. Similar expressions can be 
developed for G, G v , H , and H v . 
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2.3 Curvilinear Coordinate Transformation 


Most practical fluid flow problems of interest, are solved in domains with irregular 
shapes and boundaries. This causes difficulty in implementing the boundary conditions. 
In regions of high gradients, (shock waves, vortex regions, and shear layers), one needs 
to pack grid points in order to capture details of the flow field, and render accurate 
results. The uneven packing of grid points complicates the differencing operator. To 
avoid these difficulties, the governing equations can be transformed into a body fitted 
coordinate system, thus simplifying the numerical differencing and the implementation 
of boundary conditions. 


The curvilinear coordinate system is assumed related to the Cartesian coordinates by 
£ = £{x,y,z}, Tf = i){x,y,z}, C = C{z,y> z ) 

9,5 d_ d_ 
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The transformation matrices are given by [4]; 
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where J is the Jacobian of the transformation; 
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These formulations allow the governing equations to be written as; 
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F v , G v , and H v , are presented in Appendix A. 

2.4 Thin-Layer Navier-Stokes Equations 


(2.25) 


(2.26) 


The Navier-Stokes equations govern the motion of unsteady compressible fluid flow. 
The solution of the Navier-Stokes equations require a fine grid to capture the diffusive 
effects. Performing the computations on a fine grid requires extensive amounts of 
computer time and memory. At high Reynolds numbers, the effect of viscosity is confined 
to a thin region near solid walls where a boundary layer exists. The dominant viscous 
effects in the boundary layer arise from viscous diffusion normal to the body surface. A 
desirable approximation is to neglect the viscous terms containing derivatives in directions 
which are tangent to the body surface [26]. This assumption is often justified since the 
viscous terms containing derivatives in directions parallel to a solid boundary are usually 
substantially smaller than the terms with derivatives normal to the boundary. It would 
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also be impractical to think of a fine grid in all three directions. Viscous grids are usually 
dense along only the solid walls. Thus it makes sense to neglect the terms that are 
not properly resolved, especially if they are an order of magnitude smaller than other 

viscous terms. 

The thin-layer Navier-Stokes equations are derived from the Navier-Stokes equations 
by neglecting all cross derivatives in the viscous fluxes F v , G v , and H v . For example 
all derivatives with respect to tj and £ in the F v viscous flux are neglected. Similarly 

for G v , and H v . 

During development of the numerical algorithm, it was desired to maintain generality. 
Since it is not known a priori which direction will coincide with the solid boundary, or 
whether there will be more than one boundary surface, the thin shear layer approximation 
was applied in all three directions. The thin shear layer equations used in the developed 
algorithm are given by 
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where Q is the state vector of dependent variables. F, G, and H are the inviscid fluxes, 
described in eq. 2.25, and 
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where, <f 2 = Cj + C* + C? , and, ^ = ^(«cC*+^C»+«^C*) 

Based on the type of problem computed one or more viscous fluxes can be 

neglected [76, 96]. 


2.5 Reynolds-Averaged Navier-Stokes Equations 

Almost all flows encountered in fluid mechanics are either fully turbulent or partially 
turbulent. The nature of the flow and the purpose of the numerical study dictate the 
accuracy levels for modeling turbulent effects or the justification for neglecting turbulent 
effects completely (and simply assuming laminar flow). Turbulence enhances the rate of 
heat transfer, and alters the skin friction. Turbulence also affects the location of flow 
separation, the mechanism for separation, and the size of the separation bubble. Surface 
pressure forces, lift and drag, are also affected by the level of turbulence in the flow. 

Turbulent flows are in principle still governed by the Navier-Stokes equations, eq. 
2.23; however extremely fine grids and higher order schemes are required to resolve all 
time and length scales that accompany realistic turbulent flows. This type of computation 
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is known as Direct Numerical Simulation; it represents a challenge to today’s computers 
and numerical algorithms. Direct Numerical Simulation (DNS) has been restricted to 
low Reynolds number flows since the number of grid points required is proportional to 
the 9/4 power of Reynolds number, [97]. The significant cost of DNS calculations, even 
for simple flows, makes them impractical for current engineering applications. Perhaps 
with the development of new computer architectures, and with more parallel machines, 
DNS will become a practical approach. From that point of view, if the grid is coarse and 
must still resolve the mean details of turbulent motion, then we must resort to modeling 
the turbulent effects by superimposing them on the mean flow. At present, turbulence 
modeling forms the basis of most of the computational work in turbulent flows. 

Hinze [98] best described turbulent flow as “Turbulent fluid motion is an irregular 
condition of flow in which the various quantities show random variation with time and 
space coordinates so that statistically distinct average values can be discerned’'. Follow- 
ing the footsteps of Reynolds, we decompose the randomly changing flow properties into 
mean and fluemating components 

q = q + q (2.31) 

where q is the property being decomposed and q is the fluemating component; q is the 
mean property defined by 

to+At 

<2 ' 32) 

to 

At is a time interval which is long compared to the period of any significant turbulent 
fluctuations, but At is assumed to be short compared to the time scales associated with 
the mean flow. If we apply the decomposition procedure to all state variables in the 
Navier-Stokes equations, we get the Reynolds-Averaged Navier-Stokes equations which 
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woik well for simple incompressible flows. For compressible flows, triple correlations 
involving density appear in the equations. Favr6 [98] suggested a mass weighted 


decomposition for compressible flows to avoid the triple correlation involving density. 


The following formulations were used to decompose the flow variables in the Navier- 


Stokes equations 2.1. 



where u = u -f u, T 


a . f 

P P 

f + T, H = H + H 


note p = p + P, P = P + P 


(2.33) 


while (p + p)q = 0, 


but q / 0 

Substituting the above formulations into the Navier-Stokes equations and averaging in 
time, we get the mass averaged Navier-Stokes equation. The details of this procedure are 
presented in Appendix B. The non-dimensional mass averaged Navier-Stokes equation 


is given by 



(2.34) 


where 


F = 


pu \ 
pu 2 + p I 
< puv 
puw 

. (£ + p)u - 


0= 


p 

pu 

pv 

pw 
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A/, 


r —(h + ht){ 


«i ( 5 


(2.35) 
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2 (o du du? \ 


du i dv 

3 ^ + 37 

du i dw 

FI + 37 
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(P + PTM 


0.0 

du i dv 

3? + 37 

2 ( f\dv du 

3^ - 37 

dv i dw 

37 + 37 

f du , dv\ ,,, 2 ( ndv _ du _ dw \ , 

U 137 + 37) +”1^37 37 


c\dv du dw \ 

2ur ~ tt: --$7) 
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( dv i duA . _dT j 


H = 


i 


piU 

puw 
pvw 
pw 2 + p I 

(£ + p )™ J 


M, 


#„ = ^(P + Pt)< 

JWef 


0.0 


du i dtf 

St ?E 

37 + 37 


1(2^ - §r-Sf) 
«(fe+fe)+»(^+t)+ 


(2.38) 


where a = 


>, /!rT/ 


(2.39) 


/?(P + Pt) V^r 

p is the molecular viscosity, and pr is the eddy viscosity. The eddy viscosity is supplied 
by the turbulence model. Pr is the laminar Prandtl number and Pr T is the turbulent 
Prandtl number. For air, we take Pr = 0.72 and Prj = 0.9 

Replacing p in the Navier-Stokes equations, eq. 2.17, with p + pr and replacing 
_^-L— ^ ^ yields the modeled Reynolds-Averaged Navier-Stokes 

equations. Thus one can conclude that the mathematical formulations of the two sets 
of equation are similar. 
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2.6 Baldwin-Lomax Algebraic-Turbulence Model 

In this study, the Baldwin-Lomax turbulence model was selected to model the eddy 
viscosity. The model is a two layer eddy viscosity model which implements a simple 
algebraic expression to determine the turbulent eddy viscosity [99]. 

The inner layer eddy viscosity model is given by ; 

CB- = ^ 


where 


l = kiy 




(2.41) 



(2.42) 


and 

pw^max 

Where, )u>| is the magnitude of the vorticity; y is the normal distance from the nearest 
solid wall, Jfci is a constant equal to 0.4, A + is a constant equal to 26.0, ui max is the 
maximum vorticity in a local vorticity profile along the coordinate direction normal to 
the wall, p w is the density at the wall; and p w is the molecular viscosity at the wall. 
The original Baldwin-Lomax algebraic-turbulence model did not implement ui max in 
the formulation but rather suggested using the shear stress at the wall, t w . If there is 



Rref 

M re/ 


(2.43) 


separation, implementing the shear stress at the wall yields inaccurate values for the 
turbulence model. If there isn’t any flow separation on the wall, it can be shown that 
u m ax equals approximately t w . The second set of equations for the outer layer of the 
Baldwin-Lomax algebraic -turbulence model are given as; 


PT outer — h 2 C cp pF, wakeFkleh 


Kft / 

AW 


(2.44) 
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where 


Vmaz F i 


max 


F wake = min < 


Owk 

Ernax 


(2.45) 


with the closure constant, = 0.0168, C C p - 1.6, C w k - 1-0 for transonic flow, and 


V d if = (Vu 2 + v 2 + u> 2 ) ~ (\/« 2 + u 2 + w 2 ) ( 


(2.46) 


along the coordinate perpendicular to the surface at a particular wall location. For 
example, equation 2.46 would be applied along a constant x-surface, if x is the streamwise 
direction. The value y max corresponds to the location with F(y) = F{y) max , where 


F(y) = yM 


1 - eip {~M] 


(2.47) 


In wake regions, the exponential term for F(y) is set to zero. Alternatively, the Klebanoff 
intermittency factor is used where; 


Fkleb - 


1+5.5 


yCKLEB 


i 6 


Vmaz 


-1 


(2.48) 


and Ckleb - 0.3. 


2.7 Euler Equations 

In this section the governing equations for inviscid compressible fluid flow are 
presented. This governing set of equations is known as the Euler Equations. The Euler 
equations are considered as an approximation to the Navier-Stokes Equations. They are 
derived by neglecting the viscous and heat flux terms from the Navier-Stokes Equations, 
eq. 2.23. Euler Equations are useful models in fluid flow problems where only information 
on the pressure distribution is required. Numerical results for the Euler equations are 
important in preliminary studies and design, and the Euler equations are capable of 
capturing shock wave interactions and contact discontinuities accurately. 
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The three-dimensional, time-dependent, compressible Euler equations are written in 
general curvilinear coordinates (Z.p.C) in a non-dimensional conservative form as : 


8Q dF_ dG dH_ 
dt + d( + drj dC 


(2.49) 


where Q is the vector of dependent variables and is given by ; 


( P 



\ 


pu 

pv 


pw 

E ) 


(2.50) 


F, G, and H are the inviscid flux vectors which are functions of the state vector Q and 


are given by 


' pU ' 


' pV ' 


' pW 

puU + pit 


puV + pr)x 


puW + p(x 

< pvU + piy 

►i G - < 

pvV + prjy 

> , and H = < 

pvW + p(y 

pwU + piz 


pwV + pr) z 


pwW + p(z 

„ ( E + p)U , 


(E + p)V J 


. ( E + p)W , 


(2.51) 


The Euler equations are a set of first-order hyperbolic equations in time. For steady- 
state calculations, the mathematical nature of the equations changes from elliptic in the 
case of subsonic flow to hyperbolic for supersonic flow. 

Another set of equations that model compressible, irrotational, inviscid fluid flow 
are the potential flow equations. The potential flow equations assume that the flow is 
irrotational (y x q = 0), and do not account for vorticity and entropy changes. Thus 
they are not capable of accurately capturing curved shocks, shock interactions or contact 
discontinuities [92]. Hence the Euler equations were selected as the set of equations 
that govern the inviscid compressible fluid flow in the numerical algorithm developed 
in this study. 
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CHAPTER 3 

SPATIAL DISCRETIZATION 

A control volume approach was employed in this study. The computational domain 
was divided into a finite number of hexahedral cells. Each cell has its own volume V 
and is bounded by its surface, S. Figure 3.1 shows a typical cell in the computational 
domain. Each cell has six quadrilateral faces; each face defined by four vertices. The 


T1J 



Figure 3 3.1 Schematic of Computational Cell 


face of a hexahedral cell can collapse to an edge or even to a point, i.e., the four vertices 
may lie at the same x, y, z location. By applying the basic principles of conservation of 
mass, momentum, and energy to a stationary cell in the computational domain we extract 
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the integral form of the governing equations 


Ftl J f Qdv + I ! r ' nds = 0 

V S 


(3.1) 


where 

T = (F — F v )i + (G — G v )j + ( H — H v )k , 

(3.2) 

and n — Tixi TiyJ d - Ti^k. 

Q is the vector of conserved variables; F, F V ,G, G v , H, and H v are the flux-vectors which 
are defined in eq. 2.23; n is the outward pointing unit vector normal to the surface S, 
bounding the volume V. In the control volume approach, the state variables, Qij,k, are 
stored in the center of each cell, and are considered to be a cell average rather than a 
pointwise value at the cell center 

= 4 JJJ QdV 

K v 

The merit of the integral form of the governing equations is that it is valid everywhere in 
the computational domain, even across shocks and contact surfaces, while the differential 
form of the equations are valid only in smooth regions. The use of the conservative law 
form permits shock waves to be captured as weak solutions to the governing equations 
where the discontinuities evolve as part of the solution, and are captured within one or 



more grid cells. 


A semi-discrete form of the differential equation, eq. 2.23, can be written as 


+ (G - G v ) l j+ ^ k -{G- GtOij-i,* 

+ (H - H v ) l j k+ ^ - (H - H v ) t j k _i = 0 


(3.4) 


where A£, At?, and A ( are taken equal to unity for simplicity. If the surface integrals in 
eq. 3.1 are written as the sum of the contributions from the six faces of the computational 
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cell, then the relation between the integral form of the governing equation, eq. 3.1, and 
the semi-discrete form of the differential equation, eq. 3.4, becomes obvious. As a result, 
a geometrical interpretation of the Jacobian and metric terms of the transformation can 
be made. The Jacobian is calculated as the inverse of the cell volume; the vector y is 
the directed area of the cell interface to a / = constant coordinate direction, / = £, 77 , and 
(; 1^1 is the area of the cell interface; the unit vector consists of the direction 

cosines of the cell interface. The evaluation of the cell volume, and the directed area 
is given in Refs. [92, 100, 101]. 

The finite volume method, when coupled with an explicit scheme, has the advantage 
that the spatial discretization is decoupled from the temporal discretization. The temporal 
discretization will be discussed in detail in the next chapter, while in this chapter the 
construction of the interface flux, and the spatial discretization will be emphasized. 

The interface flux can be constructed either by a central-difference scheme or 
an upwind-scheme. Upwind schemes attempt to construct the flux by modelling the 
underlying physics, as dictated by the sign of the characteristic waves. Central-difference 
schemes compute the interface flux as an average of the two adjacent cells, disregarding 
the characteristic theory. A comparison between central-difference schemes, and upwind 
schemes, as well as how they are related, is given by Swanson and Turkel in Ref. [56]. 

The numerical procedure employed for the evaluation of the interface flux utilized a 
second-order, central-difference approximation to compute the viscous flux (F v , G v , H v ) , 
while a high resolution upwind shock capturing scheme was used to compute the inviscid 
flux ( F , G, FT). Roe’s flux -differencing scheme and Van Leer’s flux-vector splitting-scheme 
were applied to the conservative forms of the governing equations to evaluate the inviscid 
flux. 
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A Monotone Upstream-centered Scheme for Conservative Laws (MUSCL) was used 
to compute the inviscid flux at the cell interface. The state variables were extrapolated 
from the cell center to the cell interface and the interface flux was computed subsequently. 
An extrapolation of either the primitive or conservative variables can be performed. For 
the most part, the primitive variables were extrapolated in this study, because they render 
a smoother solution across shocks and slip lines and allow the use of a higher CFL 
number, as will be shown in the results section. The extrapolation operator of the state 
variables for the cell interface is based on the so-called K-scheme formulation where the 
state variables to the right and left of the interface are computed as follows 


1 = Q* + t[ — K )^« + (i + ] 

«+2 4 ( 19 ) 

Q*,1 = Qx + 1 k)V,+i + {l + k)A,+i ] 

,T 2 4 


where. A, and y, are the forward and backward differences, respectively, defined as 

At = Qt+l Qi 


( 20 ) 


V, = Qi - Qi - 1 


The value of the parameter k determines the spatial accuracy of the scheme. Table 1 


shows the conventional choices for k, the corresponding accuracy and the truncation error 


based on one-dimensional spatial analyses [101]. The parameter, <j>, is set equal to one for 
high-order extrapolations, and to zero for first-order extrapolation. Although first-order 
schemes possess good damping characteristics and allow a higher CFL number to be used, 
they are too diffusive. In the present study, k = 1/3, k = -1, and k = 0 have been used. 


High-order-accurate upwind schemes produce spurious oscillations and may develop 
instabilities near shocks and contact discontinuities. These oscillations can be reduced by 
the use of some kind of limiting procedures called limiters. This is achieved by imposing 
a constraint on the gradients of the dependent variables or on the flux functions. Several 
limiters are available in the literature [72, 102-104]. In the present study the min-mod 
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limiter and the differentiable Vanalbda limiter have been employed. The two limiters are 
discussed in Appendix C. The limiters prevent numerical oscillations, but at the same 
time reduce the accuracy of the scheme to first order near discontinuities. The main 
disadvantage of limiters is that they lead to a limit cycle, with no apparent convergence. 
The residual will converge to a certain level and then “hang up”. To overcome this 
problem, the second-order strictly upwind formulation, k = -1, was used which does 
not require a limiter. 

In summary, the present algorithm has employed a cell-centered, finite-volume ap- 
proach with the state variables at the cell interface determined by the MUSCL interpo- 
lation. The viscous and diffusion terms, at the cell interface, were discretized using a 
second-order, accurate, central-differencing scheme, while a higher-order-upwind-scheme 
(Roe’s flux-differencing scheme [66] or van Leer’s flux-vector splitting- scheme [59]) was 
used to construct the inviscid flux. Both Roe’s flux-differencing scheme and van Leer’s 
flux-vector splitting-scheme are capable of capturing relatively strong stationary —shocks 
within one or two interior cells— if the shock is reasonably aligned with the grid. It 
should be emphasized that the present algorithm employed a sequence of one-dimensional 
operators in all three coordinate directions £, rj, and £■ Thus the use of highly skewed, 
non-orthogonal grids should be avoided, if possible, because the one-dimensional opera- 


Table 1 Values of k and the Corresponding Truncation Error. 
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tors assume that the waves interact normal to cell interfaces. This method of extending 
one-dimensional schemes to multi-dimensions was found to be quite satisfactory for the 
test cases investigated in this study. A truly multi-dimensional approach is still in its 
infancy, and is computationally expensive [5, 105, 106]. 

3.1 van Leer’s Flux-Vector Splitting 

A prevalent way to introduce upwinding into the governing systems of hyperbolic 
conservation laws is to split the flux according to the characteristic speed (q, q±c)). In 
this study, van Leer’s flux-vector splitting has been employed, because it yields sharp, 
crisp shock surface. A disadvantage of van Leer’s scheme is that it smears slip lines 
because it ignores the linear waves (entropy and shear waves) [107]. A brief summary 
of the scheme is presented in this section. For more information about the scheme, the 
interested reader should review the original paper of van Leer [59]. 

Following Ref. [59], the flux-vectors F, G, and H can each be split into two vectors, 
a forward flux-vector, based on non-negative eigenvalues, and a backward flux-vector 
based on non-positive eigenvalues. 

F = F + + F“, G =G + + G", H = H + + H~ (3.7) 


For local supersonic Mach numbers: 

Mi > 1.0, Ff = F h Ff — 0 


(3.8) 


Mi < -1.0, F, + = 0. F,- = F t 

where / = £, rj, and ( to indicate the three coordinate directions. 

For subsonic local Mach numbers, |A//| < 1.0 (in general notation for body fitted 
coordinates [53]), a local scaled contravariant velocity component u\ is defined as 


IjU + lyV + IzW 

yjf- + Ty + /“ 




(3.9) 
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where the local Mach number is given as 


ui 

Mi — — 

a 


and a is the local speed of sound. The fluxes are: 


where, 




i/* 

j J ma 


i x {-ui i 2a)/7 + u 

ly{—ui ± 2 a)j~f + v 

l z (-ui ± 2a)/7 + w 

ft 


F* 

t energy 


ln = 


In 


\J~ii + l l + 


n = x, y, and z 


/mass = ±pa-{Ml±l)‘ 


and 


r± _ r -3fi?±2gtiia+2fl 2 u 2 +t» 2 +w 2 
J energy [ -yi — 1 ' 2 


Here; 


F{ = F F n = G F C = H 


U£ = U U v = V Uq — W 


and 


0 = 7 - 1 

The “+” indicates a forward flux and the indicates a backward flux. 


( 3 . 10 ) 


( 3 . 11 ) 


( 3 . 12 ) 


( 3 . 13 ) 


( 3 . 14 ) 


( 3 . 15 ) 


( 3 . 17 ) 
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Applying the split flux to the semi-discrete form of the governing equations, eq 2.27, 


gives: 


dQi,j,k 

dt 


+ (F + + F- - F v ) i+ , jJc - (F* + F-- F„),_ i M 
+ (G+ + G- - G,) l j+ , t - (G+ + G-- G v ) t j _i k 
+ + H~ - W„) lji+ , - (ff + + H~ - j 


0 


(3.18) 


The present formulation, when applied to transonic and low supersonic flows, does 
not require the use of flux limiters for essentially oscillation free shocks. This was pointed 
out by Anderson, Thomas, and van Leer [58], von Lavante and Haertl [108], Melson and 
von Lavante [109], and Cannizzaro, von Lavante, and Melson [110] and was explained 
in more detail by van Leer [59] 

3.2 Roe’s Flux-Difference Splitting 

Roe’s flux-difference splitting is an upwind scheme that approximates the Riemann 
problem at an interface between two cells by Roe’s averaging procedure [66]. Roe’s 
scheme provides an exact solution to an approximate Riemann problem, and is capable 
of handling slip lines with less smearing. The idea of advancing the solution to the next 
time-level by solving a Riemann problem at each cell interface was first introduced by 
Gudonov [64], The Riemann problem and the different waves associated with it are 
illustrated in Fig. (3.2). A good review of the different Riemann solvers is given in 
Refs. [6, 72, 111]. 
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Figure 3.2 The Riemann Problem. 

Roe’s scheme is used widely in practical applications of computational fluid dynamics 
because it results in an efficient and accurate computation form. A comparison of 
different numerical flux formulas for the Euler and Navier-Stokes equations. Ref. [107], 
recommended the use of the Roe scheme. That recommendation was based on the fact that 
the interface flux computed by the Roe scheme includes information about all different 
waves — linear and non-linear — by which the neighboring cells interact. The scheme 
gives good results when shocks, contact discontinuities, and slip lines are aligned with 
the grid. The main advantage of the scheme is that it returns the exact solution whenever 
qL and £>* tie on opposite sides of shocks and contact discontinuities. However, Roe’s 
scheme can also represent an expansion shock due to the lack of a naturally constructed 
entropy condition. 

The interface flux-vector F/ is evaluated as 


(f, + i) ( = - 1 - <?£.)] (3.i9) 


Fi and Fr are the flux-vectors computed from the left and right states, and A is the 
Roe averaged flux Jacobian matrix 


A = 



(3.20) 
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where 



(3.21) 


and 


= S { |A|S s - 


-1 


(3.22) 


and 5 ’” 1 are the right and left eigenvectors of Roe’s averaged Jacobian matrix A; A 
is the diagonal matrix of the eigenvalues of A. The Roe averaged state is computed as 


P = yfpRPL 


u = 


v — 


w = 


H = 


a = 


ULy/PL_ + UR^/pR_ 
yfPL + y/PR 
V_Ly/PL_ + VRy/PR 
y/PL + y/PR 
WL\/PL_ + ^Ry/PR 
yfpL + y/PR 

thyfpl^ Hr^/pr_ 

y/PL + y/PR 


'(7-1) 


H — ^{w 2 + £’ 2 + ti> 2 } 


(3.23) 


The tildes refer to the Roe averaged quantities. The last term in eq. 3.19, i.e. 
|A|(Qfl - is a damping term due to the upwind character of the scheme and is 
given in detail in Ref. [55] as 


A 


( Qr - Ql) = 


0 4 

UO4 + ^05 + 06 
VQ4 + lyCtS + Orj 
lZ ’04 + / z O 5 + 08 
a 5 


(3.24) 


where. 


a 5 = 


V 04 + u /05 + uo 6 + voti + u>o 8 — 


~ 9 

a\a 
7 - 1 


(3.25) 
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and 
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\ui\p(Aw - l z Aui) 


(3.26) 


with Vf = J l l + l y + l l . for / = f, r?, or C- The A refers to the difference between 
the state variables on the left and right sides of the cell interface, such as A p = p fl - p L 
Here; 

F c = H (3.27) 


Fc = F 


F n = G 


and 


uj = u = v U C = w (3.28) 

For more information about the Roe scheme, the interested reader should refer to the 
original work by Roe, Refs. [65, 66]. 
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CHAPTER 4 
TIME INTEGRATION 

In this study emphasis was placed on the simulation of complex three-dimensional 
steady, and unsteady fluid flow problems. The selection of a particular type of time 
integration technique, whether implicit or explicit, will determine the characteristics of 
the numerical method used to investigate the fluid problem of concern. In this study, 
explicit time-stepping schemes were used to construct the algorithm for solving the time- 
dependent, compressible, Euler, and Navier-Stokes equations. 

There are a large number of explicit schemes that have either been used previously or 
are still in use for solving the compressible flow equations. The desired numerical method 
should be simple, robust, have effective damping of high frequency errors (necessary 
for multigrid), have low dispersion (low phase errors will reduce spurious oscillations 
and result in faster convergence rates), low levels of numerical dissipation for accurate 
predictions of viscous effects and it should maintain high resolution on stretched grids. 
Programing simplicity is another important issue, since the goal is to implement the 
time-stepping scheme in a multi-block code, and on massively parallel machines. 

However, no generic time-stepping scheme was found that satisfied all of the re- 
quirements, thus it was decided to develop explicit time-stepping schemes that could be 
tailored to our needs. Two similar but distinct time-stepping schemes were developed 
for the purpose of solving the compressible, time-dependent, governing set of equations. 
The two schemes were the multistage time-stepping scheme, and the Predictor-Corrector 
Scheme. Both schemes are explicit, but they are distinct since each of the techniques 
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utilizes a different operator to compute the flux at cell interfaces. In order to understand, 
compare, and assess the two time-stepping techniques, they were applied initially to the 
model wave equation. There a conventional Fourier stability analysis could be carried 
out, yielding amplification factors, and the stability characteristics for each scheme. Dis- 
cussion of the evaluation of the two schemes will be presented and discussed in this 
chapter. 


4.1 Multistage Time-Stepping Scheme. 

Modified Runge-Kutta methods, with standard coefficients, have been rather success- 
ful when used in combination with central-difference, spatial-discretization techniques. 
Unfortunately, they perform very poorly with upwind differencing schemes. In this sec- 
tion an attempt has been made to modify the standard coefficients to achieve better per- 
formance, resulting simultaneously in schemes that are, in general, of reduced accuracy 
in time. To explore the damping properties and extend the stability limits of the explicit 
multistage scheme, the scheme was first applied to the ordinary differential equation 


dq 

It ~ ~' q 


and which has the analytical solution: 


w 


here z > 0. 


(4.1) 


q = K 


(4.2) 


Here, q 0 is the initial value of q at t = t„. 

A Taylor series expansion for q n+l around (f gives 


cit dt • 


9 = V T — “ ' T oi ' dt 3 3! 

Substituting for the derivatives of q in the above equation 


(4.3) 


<7 n+1 = q U 


l + (-:Af) + (-:Af)‘ + (-:Af) + 


(4.4) 
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The above equation is the foundation of a numerical integration scheme known as the 
modified Runge-Kutta Scheme [27]. 

Consider now the model wave equation in the form 


dq dq 

S + “& =0 

where a is the wave speed which is assumed to be real and positive. Then 


(4.5) 


dq dq 

-T- = -a-T- = -d(q) 

at ox 


(4.6) 


where R(q) represents the right hand side of eq. 4.6. If we assume that q can be 
represented in an exponential form, then 


q 


= <?o e 


-z(t-t a ) 


(4.7) 


where q 0 is the initial value of q at t = t 0 . Therefore, 


- -a , e -*V-*o) 

dt ~ q ° 


(4.8) 


8q 

dt 


-zq 


(4.9) 


The multistage explicit time-stepping can be used to advance eq. 4.9 in time from time 
step n to n+1 in the following way 

q = q 

q ] = q" + 0] A tR(q°) 


q- = q n +o 2 A tR(q l ) 
q™ = q n + a m AtR{q m - 1 ) 


(4.10) 


q 


u+1 


= q m 


41 



Comparing eq. 4.9 and eq. 4.6 yields 


/!(,) = (4.11) 

OX 

If we define 

P = At x r, (4-12) 


combine P with eq. 4.11 and substitute into eq. 4.10, we get 

q n+l =q n (\ + a m P + o m o m _,P 2 + -. + a, 02 ....a m P m ) (4.13) 

By comparing the terms in eq. 4.13 and eq. 4.4, we can determine the temporal accuracy 
of the multistage scheme. The scheme will be first order in time if 

a m = l. (4-14) 


The scheme will be second-order in time if 


a ni = 1. and o m -i = 


(4.15) 


The scheme will be third-order in time if 


a m = 1. o,„_i — — . and <*m-2 


3' 


(4.16) 


The scheme will be fourth-order in time if 


a m = 1. Om-l = “ ■ o ,„_2 — - and. a m -3 — ^ 


(4.17) 


and so on. 

One can continue this progression and arrive at higher order schemes. It should be 
emphasized that the leading coefficient o m should always be 1.0 for the scheme to be 
at least first order in time. 
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The amplification factor of the multistage explicit time-stepping scheme, G, can be 
derived from eq. 4.13 and will take the following form 

C — — = 1 4- ft]/ 5 d-aio^f’’ + ••••"(' oi (4. 18) 

q n 

The stability and damping properties of the multistage scheme are associated with 
complex polynomial G. G is a function of the coefficients a s and of P. The complex 
function P is a function of the spatial operator used to interpolate q at the cell interface. 
Thus the stability properties of the multistage scheme are tied to the spatial operator used 
to compute the flux at the cell interface. 

In this study a control volume approach was implemented where the spatial dis- 


cretization of the wave equation takes the form: 


da 

R {q) = a - = a 


*+* " q '-h 


Ax 


(4.19) 


The extrapolation of the state variables to the cell interface is base on the so-called 
K-scheme where. 


A " , = q" 1 + -[(1 — *)A, + (1 + «)Vi] 

2 4 

(/ n _, = Cl 1 + t[( 1 - *)A,-1 + (1 + K)Vt-i] 

l 2 4 


(4.20) 


such that 


— q i qi — i and \y j i— i ^/i 


(4.21) 


As mentioned in the previous chapter k determines the spatial accuracy of the scheme; 
K = -1 is a fully upwind second-order accurate scheme; k = 0 is an upwind biased, 
second-order Fromm scheme; h = 1/3 is an upwind, biased third-order accurate scheme, 
and k = 1 is a second-order accurate, central-difference scheme. The first-order scheme 
is obtained by setting / to zero. For simplicity, in the present stability analysis the limiter 
was not included in the K-scheme. 
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If we now assume the data to be harmonic 


< 7 j — tfx=j Ax — Qo£ 


13 ] 


(4.22) 


where ,3 is the spatial wave number ranging from 0 to tt, and / is \J— 1. Values of 3 
between tt/ 2 and n are considered to be high frequencies. Combining eq. 4.19, eq. 4.7 
with eq. 4.20 yields 


P = 




’ 1 - K ( 

1 - t ~ l& \ 

4 ' 

1 e j 



(4.23) 


Fourier transform of the spatial-operator (P) is a function of the CFL number and the 
wave number, P {CFL, 3). The expression for the amplification factor, given by eq. 4. 1 8, 
defines the stability region of the scheme. The stability of the multistage scheme requires 
that the modulus of the amplification factor IGI be less than unity. This expression gets 
complicated if we attempt to substitute the expression for the complex polynomials, P, 
into the expression for G and define the stability region of the scheme analytically. An 
alternative way to determine the stability region of the scheme is to plot the modulus 
of the amplification factor for the multistage explicit scheme and identify the stability 


limit graphically. 

The stability of the multistage scheme depends on the complex polynomial P and 
the coefficients a s of the multistage scheme. The locus of the Fourier transform, P, 
superimposed on the contours of the amplification factor can be used to optimize the 
coefficients of the explicit, multistage time-stepping scheme to better suit the upwind 
schemes, and achieve better rates of convergence to steady-state. The modulus of 
amplification factor \C\ with the locus of Fourier transform for a first-order, four— stage 
standard Runge-Kutta scheme are shown in Fig. 4.1 The influence of the coefficients 
on the contours of the amplification factor had to be fully understood to facilitate the 
selection of the optimum coefficient set. Optimization of the coefficients, a s , was earned 
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Figure 4.1 Contours of Modulus of Amplification Factor, |G|, and Locus of Fourier 
Transform, P for a First-Order, Four-Stage Runge-Kutta Scheme, CFL = 2.0 


out by displaying the stability plots on a computer terminal. The changes in the shape 
of the contours of IGI were observed in real time as the coefficients were changed. The 
“islands” of low values of IGI correspond to the roots of eq. 4.18. The main purpose 
of the optimization was to find a combination of the coefficients, ot s , such that, for the 
largest possible CFL values, there would be good high frequency damping (low values 
of IGI) over a large range of CFL. That is, the optimal coefficients should maximize the 
size of the islands, and make them as close to the real axis as possible. The optimization 
was performed for the two-stage, three-stage, and four-stage schemes. For each of the 
mentioned schemes the optimization was conducted for four different spatial operators: 
first-order; second-order fully upwind (k =-l); second-order Fromm Scheme (k = 0); 
and third-order upwind biased (k = 1/3). Tables 1—4 list the optimized coefficients for 
the spatial operators mentioned above. 
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Table 2 Multistage Coefficients for First-Order Scheme 


Number of Stages 

Multistage Coefficients 



01 

a 2 

03 

C*4 

Two-Stage Scheme 

1.0 

0.22 



Three-Stage Scheme 

1.0 

0.325 

0.105 


Four-Stage Scheme 

1.0 

0.34 

0.152 

0.056 


Table 3 Multistage Coefficients for Second-Order Fully Upwind Scheme 


Number of Stages 

Multistage Coefficients 



C*1 

« 2 

q 3 

«4 

Two-Stage Scheme 

1.0 

0.22 



Three- Stage Scheme 

1.0 

0.4 

0.15 


Four-Stage Scheme 

1.0 

0.42 

0.24 

0.091 


Table 4 Multistage Coefficients for Second-Order Fromm Scheme 


Number of Stages 

Multistage Coefficients 



ai 

a 2 

a 3 

04 

Two-Stage Scheme 

1.0 

0.42 



Three-Stage Scheme 

1.0 

0.44 

0.21 


Four-Stage Scheme 

1.0 

0.46 

0.255 

0.11 


Table 5 Multistage Coefficients for Third-Order Upwind Biased Scheme 


Number of Stages 

Multistage Coefficients 



ttj 

a 2 

<>3 

04 

Two-Stage Scheme 

1.0 

0.46 



Three-Stage Scheme 

1.0 

0.48 

0.22 


Four-Stage Scheme 

1.0 

0.44 

0.26 

0.135 
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The simplest schemes to optimize were the two-stage versions, since only one 
coefficient can be selected freely. The first coefficient is always equal to unity to 
ensure that the scheme is at least first-order accurate in time. The challenge to optimize 
the coefficients of the explicit multistage scheme increased by increasing the number of 
stages, since the number of coefficients to be optimized increased. The most challenging 
scheme to optimize was the four-stage scheme since the optimum combination of three 
coefficients has to be found. The modulus of the amplification factor IGI with the locus 
of its Fourier transform (of the spatial operator, P, corresponding to a maximum CFL) 
number is presented in Figs. 4.2, 4.3, and 4.4, for all the spatial operators used in this 
study. The resulting stability plots will be shown only in the second quadrant (upper 
half of the negative real part of the complex polynomial P) since they are symmetric 
with respect to the real axis. Figures 4.5, 4.6, and 4.7 represent the magnitude of the 
modulus of the amplification factor IGI as a function of the spatial wave number, 8, and 
the CFL number. By displaying the two sets of plots for a particular multistage scheme, 
the stability region and the damping properties of the scheme can be fully displayed. 
By increasing the number of stages, we are able to increase the CFL number, as shown 
in Figs. 4.5, 4.6, and 4.7. The time to perform a four-stage explicit scheme is twice 
that for a two-stage scheme. On the other hand the CFL number increased form 2 to 
4, comparing with the case of a first-order scheme. This conclusion is also valid for 
the remaining spatial operators, as shown in the stability plots. The main advantage of 
going to a higher number of stages was the good damping characteristics for high wave 
numbers. Considering the results of the stability analysis, the most promising schemes 
of practical importance were the Four-Stage Fromm Scheme (k = 0) and the third-order 
upwind biased Scheme (k. = 1/3). The Fromm Scheme was preferred due to its low 
numerical dispersion, demonstrated by results with the least oscillations around shocks. 
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It is important to point out that, when multigrid acceleration techniques are imple- 
mented the desire for maximum CFL number is not as important as the high frequency 
damping requirement. The high frequency damping (or lack of it) will affect the rate of 
convergence to steady-state more significantly than the CFL number. The choice of the 
optimum CFL number, when utilizing multigrid acceleration techniques, should be based 
on how well high frequencies are damped. 

It should be mentioned here that, in a parallel effort, van Leer, Tai, and Powell [35], 
and Gaffney [112], also tried to optimize the Runge-Kutta coefficients for applications 
with the upwind methods. The van Leer, Tai, and Powell approach was somewhat 
different than the work discussed previously. Their approach assumed that a genuine 
and practical multi-dimensional characteristic formulation of the Euler equations could 
be found, and then they optimized the Runge-Kutta coefficients for only one value of the 
CFL number. They argued that each wave would propagate at its optimum CFL-number. 
Unfortunately, there is no such formulation for three dimensional cases. Generally, the 
maximum CFL numbers, for the van Leer, Tai, and Powell approach, were lower, and 
the damping was effective over a narrower range of CFL numbers. 


48 




Re(P) 

(a) First-Order, CFL = 2.0 



Re(P) 
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Figure 4.2 Contours of Modulus of Amplification Factor, |G| , and Locus of Fourier 


Transform of Spatial-Operator, P for Two-Stage Schemes. (Continued 
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(c) Fromm, Second-Order, CFL =1.0 
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(d) Third-Order, CFL = 1.0 


Figure 4.2 Contours of Modulus of Amplification Factor, |G| , and Locus of 
Fourier Transform of Spatial-Operator, P for Two-Stage Schemes. 
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Figure 4.3 Contours of Modulus of Amplification Factor, |G|, and Locus of 
Fourier Transform of Spatial-Operator, P, , for Three-Stage Schemes. 
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(b) Second-Order, CFL = 1.25 

Figure 4.4 Contours of Modulus of Amplification Factor, |G|, and Locus of Fourier 
Transform of Spatial-Operator, P, for Four-Stage Schemes. (Continued . . . ) 
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(c) Fromm, Second-Order, CFL = 2.1 
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(d) Third-Order, CFL = 2.0 


Figure 4.4 Contours of Modulus of Amplification Factor, |G'|, and Locus of 
Fourier Transform of Spatial-Operator, P, , for Four-Stage Schemes. 
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(a) First-Order 
1 

IGI 


(b) Second-Order 

Figure 4.5 Modulus of Amplification Factor, |G'|, as a Function of the Spatial Wave 
Number, 3 , and the CFL number, for Two-Stage Schemes. (Continued . . . ) 
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Figure 4.5 Modulus of Amplification Factor, |G|, as a Function of the Spatial 
Wave Number, J, and the CFL number, for Two-Stage Schemes. 
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(a) First-Order 



(b) Second-Order 

Figure 4.6 Modulus of Amplification Factor, |G’|, as a Function of the Spatial Wave 
Number, 8 , and the CFL number, for Three-Stage Schemes. (Continued . . . ) 
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(c) Fromm, Second-Order 



(d) Third-Order 


Figure 4.6 Modulus of Amplification Factor, |G|, as a Function of the Spatial 
Wave Number, $, and the CFL number, for Three-Stage Schemes. 
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Id 


(c) Fromm, Second-Order 

Figure 4.7 Modulus of Amplification Factor, |G|, as a Function of the Spatial Wave 
Number, 8 , and the CFL number, for Four-Stage Schemes. (Continued . . . ) 
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Figure 4.7 Modulus of Amplification Factor, |G’|, as a Function of the Spatial 
Wave Number, .7, and the CFL number, for Four-Stage Schemes. 
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4.2 Predictor-Corrector Schemes. 


The main goal was to develop a two-stage explicit scheme which had good damping 
qualities and could allow the use of large CFL numbers. In the previous section we 
presented two-stage schemes for different spatial operators. These schemes offered good 
damping qualities but the maximum allowable CFL number was less than unity, except 
for the first-order spatial discretization , as shown in Figs. 4.5, 4.6, and 4.7. It was decided 
then to construct a two-stage scheme where the first stage (predictor step) required the use 
of a first-order spatial operator and the second stage (corrector step) utilized a second- 
order spatial operator. The scheme is a modified version of the upwind scheme of 
Warming and Beam [18], and is named 1-2 scheme. The two steps are given by: 
Predictor Step: 


4 * 


7,-i 

2 


= 7 , 


7,-1 


q] = qf - A t * R(q ») 


(4.24) 


Corrector Step: 


( , r l+ i — 7" +9(9/ — 7 ,"-i ) 

. = 7,-1 + ^(7,-1 - 7, ”- 2 ) 


(4.25) 


q» + ' = qf - At * R(q 2 ) 
where, R(q) is the residual and is given in eq. 4.19 

A stability analysis similar to the one in the previous section, was performed. The 
only difference was that the type of the spatial operator used for the extrapolation to 
the cell interface was different. The resulting plot of the magnitude of the spatial wave 
number, /?, between 0 and 7r and the CFL number between 0 and 2 is shown in Fig. 4.8. 
The figure shows clearly that the stability limit is higher than for the two-stage schemes 
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introduced in the previous section. The scheme is still second-order accurate in space and 
satisfies the “shift condition”. Introducing the first-order spatial operator in the predictor 
step allowed us to increase the CFL number to a value of 2. The scheme is a very simple 
and efficient method consisting of only two steps. 

There were, however, some problems with the scheme. Since the two steps are 
different, the steady-state result, depended on the time step. This phenomena was 
observed in only a few cases, represented by convergence to a residual that was larger 
than “machine zero”. The second problem with this scheme was the limiters. None of 
the flux limiters tested in this scheme converged more than two orders of magnitude. 
The reason most probably is due to the mixing of time levels in the extrapolation of the 
variables to the cell interface in the corrector step. Never the less the resulting flow fields 
agreed well with other, fully converged numerical results. 

A more serious problem is the increase of the damping factor to 1.0 at high frequencies 
for CFL number 1, as shown in Fig. 4.8. In the scalar case, the CFL number can be kept 
at its optimum value of 1.7. But in the case of the Euler or the Navier-Stokes equations 
there are three distinct eigenvalues in each direction. Typically, local time-stepping can 
be implemented, where each cell is advanced at its optimum time step. The time step is 
a function of the corresponding eigenvalues at that cell. One or more eigenvalues might 
correspond to the CFL number ranges with minimum or no damping. This was manifested 
by the lack of convergence of the 1-2 Scheme when utilized in a multigrid procedure. 
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For multigrid applications, this scheme was modified by making the predictor step 
second-order accurate. Thus, the predictor step equation 4.24 is replaced by the fol- 
lowing: 


( 7, 1 + i - <?" + 2^«" - 9»”-i) 
<?'_! = <71-1 + ^ (<7t-l _ <7 i"-2) 
q}=q!-A t*R(q l ) 


(4.26) 


The corrector step is identical to the corrector step in the 1-2 Scheme and will only be 
repeated here for the sake of completeness. 

Corrector Step: 


(4.27) 


</ ] + i - <7 r + 77 (< 7 i ~ ( ii- 1 ) 

<7, 2 _i = <7"-i + 2 ^' 1_1 ~ 

< 7 r +I = 9 ” - At * 

The plot of the damping characteristics of this 2-2 Scheme is shown in Fig. 4.8. The 
maximum stable CFL number is now only 1.0, but its high frequency damping is improved 
significantly. The 2-2 Scheme performed better with multigrid. 


One of the main advantages of using the explicit time-stepping schemes was their 
simplicity. They could be extended to higher dimensions easily. They fit into a multi- 
block environment naturally. They were easily implemented on massively parallel 
machines, such as CM-2. The main drawback was the restriction on the time step, 
especially on highly stretched grids (viscous grids). However, by utilizing acceleration 
techniques, local time-stepping, implicit residual smoothing and multigrid techniques it 
was possible to overcome this drawback and increase the stability region of the scheme, 
as will be shown in the next chapter. 
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IGI 


(a) 1 -2-Predictor-Corrector Scheme 




(b) 2-2-Predictor-Corrector Scheme 

Figure 4.8 Modulus of Amplification Factor, |G|, as a Function of the Spatial 
Wave Number, 3 , and the CFL Number, for Predictor-Corrector Schemes. 
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CHAPTER 5 

ACCELERATING TECHNIQUES 


Local time -stepping. Implicit Residual Smoothing, and Full Approximation Storage 
(FAS) multigrid procedure, were employed in the present algorithm to remove the stiffness 
from the governing set of equations and to accelerate the rate of convergence to the 
steady-state solution. These techniques were integrated with the explicit-multi-stage time- 
stepping scheme, discussed in the previous chapter, to enhance the overall computing 
efficiency and performance of the algorithm. For unsteady flows, the accelerating 
techniques employed cannot be applied and a global time step was used. 


5.1 Local Time-Stepping 


Local time-stepping allows each cell to advance in time by the maximum allowable 
local time step, as dictated by the local stability requirements. This process allows 
faster signal propagation through the computation domain, relaxes the stiffness of the 
governing equations, and hence increases the rate of convergence to steady-state. An 
accurate estimation of the allowable time step is of paramount importance if a robust 
algorithm is desired. The time step used in this study was based on both a convection 
and a diffusion stability limit, and is given by: 


where 


At 

CFL 

A ti 

CFL 


1 1 1 J_ J_ 1 1 

^ A + A t v + At<; + A t V £ + Atjj A 


( 5 . 1 ) 


At? 




>>?,= 


U l 

t 

Pi 


+ a yj 0 J 
<pl[max 


j, -pA } + jdWs.1 + iw + iv»i) 


<t>i = l 2 x + % + /J; / = €,»?. and C 
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The first three terms on the right hand side of equation 5.1 result from stability limitations 
on the inviscid flux, while the last three terms are due to stability limitations on the 
diffusion terms. The expression for the time step is based on the stability analysis study 
presented in [33] . The diffusion limit on the time step » makes the scheme more 
robust on fine grids, and in boundary-layer type flows [113]. 

For unsteady flows a global time step was used which was required to be the smallest 
maximum time step calculated within the computational domain. 

5.2 Implicit Residual Smoothing 

Implicit residual smoothing extends the stability limit, and improves the damping 
properties of the multistage time-stepping scheme. Lerat [38] introduced the idea of 
residual smoothing for the Lax-Wendroff scheme. Jameson and Baker [29] applied 
the idea of an implicit residual smoothing in conjunction with modified Runge-Kutta 
schemes. This procedure has been developed further by [28, 33, 40-42], by employing 
a central-implicit-residual-smoothing operator. The use of an upwind-residual-smoothing 
operator was employed also [43, 44], The smoothing operator modified the basic k-stage 
explicit scheme of eq. 4.10 in the following manner 

dQ d(F-F v ) d(G-G v ) d(H - H v )' 

dt ~ [ d£ drj + d( . 

Q° = Q n 

Q 1 = Q a -a l R'{R{Q 0 )) 

Q 2 = Q n -a 2 R*{R{Q 1 )) (5 ' 2) 

Q k = Q n -a k R* ( R(Q k ~ ')) 

Q n+1 = Q k 
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The smoothed residual R* is a function of the unsmoothed residual R. The smoothing 
operator can be either an explicit or an implicit operator. In the present study an implicit 
operator was employed. For three-dimensional flow, the smoothing operator is given 
by: 

(1 - e*%)( 1 - t v d nv ){ 1 - e c d c< )R* = R (5.3) 

where e^, tq, and are the smoothing coefficients for the £, 77, and £ coordinate 
directions respectively. The second-order-central difference operators are 6gg, bqq, and 
6^; where 

d^R = Ri+\j,k ~ ,j,k + R-i—i, j,k ( 5 - 4 ) 

Similar expressions for Sqq, and apply. An efficient tri-diagonal solver, the Thomas 
Algorithm, was applied sequentially in all three directions to evaluate the smoothed 
residual R * such that: 

(l — e{S{{)R^ - - R 

(1 = R* (5.5) 

(1 -e c 6 cc )/?‘ = R" 

Following the same guidelines and notation used in performing the stability analysis 
in chapter four, an analytical investigation of the one-dimensional wave equation was 
conducted to study the effect of the residual smoothing operator on the stability limit 
of the basic explicit time-stepping scheme. The implicit residual smoothing operator, eq 
5.3, for a one-dimensional problem, in the x-direction, is given by 

(1 -t x 6 xx )R* = R (5.6) 
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which yields the following form when expanded 


— txR*—\ + (1 + 2 t x )R t — f-xRi — i — R i 


(5.7) 


If we assume that R* = R*(t)e lu,z and R = R{t)e Iu,x and substitute in eq. 5.7, we get 


R* = 


R 


1 + 2t x {\ — cos((3)} 


where /? = ujx 


(5.8) 


The above equation shows that the smoothed residual is a function of the unsmoothed 
residual, the smoothing coefficient, and the spatial wave number. 

The implicit residual smoothing performed well, when combined with the modified 
Runge-Kutta scheme, and the 2-2 scheme. The implicit residual operator damped the 
high frequency errors and allowed the use of a higher CFL number which improved 
the rate of convergence to steady-state. On the other hand, the overall performance 
was rather disappointing when combined with the 1-2 Scheme, and no gain from using 
the implicit operator was achieved. The plot of the damping characteristics for the 2- 
2 Scheme, when combined with the smoothing operator and a smoothing coefficient of 
0.5 is shown in Fig. 5.1. Comparing Fig. 5.1 and Fig. 4.8 for the unsmoothed scheme 
clearly demonstrates that residual smoothing not only increased the stability limit of the 
scheme and allowed the use of a higher CFL number, but it also provided good high 
frequency damping. 

In the early stages of the present work, a constant scalar residual smoothing coefficient 
was used. The smoothing operator was employed after every stage of the time-stepping 
scheme and was activated uniformly in all three directions. The value of the residual 
smoothing was selected to be between 0.1 and 0.5, depending on the case investigated, the 
computational grid, and the CFL number. A higher value for the smoothing coefficient 
was used in a coordinate direction where the grid was highly stretched. Increasing the 
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value of the smoothing coefficient allowed the use of a higher CFL number. Increasing 
the CFL number by a factor of two, usually gave the best rate of convergence. It should 
be emphasized that changing the value of the smoothing coefficients changed the shape 
of the amplification map as well as the stability range of the scheme. 


In the next stage of development, an adaptive implicit residual smoothing technique 
was employed. This procedure was originally suggested by Martinelli [40], and de- 
veloped further in by Swanson, Turkel, and White [33] for two-dimensional, central- 
difference schemes. It was extended subsequently to three-dimensions, and yields the 
following expression for the smoothing coefficients [42, 114]: 



Figure 5.1 Modulus of Amplification Factor as a Function of the Spatial Wave Number 
and the CFL for the Predictor-Corrector Scheme with Residual Smoother (e = 0.5). 
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(7FI (A^ + + A c )) 

( CFL* A„ 

{ CFL (Xr, + Uh + \)) 


{ CFU_ \ 

V CFL (A c + jk(A^ + A,)) 





(5.9) 


Here, e^, tq, and are the adaptive residual smoothing coefficients which are 

CFL*' 

functions of the grid aspect ratio and the spectral radii A £,\ v ,and A^; Wr is tht 
ratio of the CFL number of the smoothed scheme to that of the basic explicit scheme. 
Increasing the ratio of C pfp L [ >2.0 caused the high frequency damping of the scheme to 
vanish, which was detrimental to multigrid convergence. The smoothing operator was 
applied after every stage of the Runge-Kutta time-stepping scheme. 

5.3 Multigrid Method 

The multigrid acceleration technique has been employed in the present work to 
augment the time-stepping schemes, discussed in chapter four, and to enhance the 
performance of the developed algorithm. Multigrid is still in its infancy, and a great 
deal remains to be learned about its performance. The multigrid acceleration technique 
was developed originally by Fedorenko [45, 46] in 1961. It was further developed 
by Brandt [47] and applied to an elliptic set of equations. The work by Brandt and 
many others has led to the popular use of multigrid by many in the fields of applied 
mathematics and computational engineering. The basis for multigrid is the use of 
successively coarser grids to calculate corrections to the solution of a partial differential 
equation (or set of partial differential equations) on a ‘fine’ mesh. These corrections 
reduce the low frequency components of the error in the fine-grid solution. Since the 
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coarse-grids contain significantly fewer points than the fine-grid, less work is required to 
perform a computation there than on the fine-grid. Excellent discussions concerning the 
development of the multigrid technique can be found elsewhere [48-50]. The multigrid 
has been used successfully for solving the potential, Euler, and Navier-Stokes equations 
[51-53, 55] 

The development and implementation of multigrid for linear problems is described by 
Briggs [115]. Unfortunately, many problems in engineering are described by non-linear 
equations or sets of equations. This is particularly true for computational fluid dynamics. 
Because of the non-linear nature of the equations, the Full Approximation Storage (FAS) 
multigrid procedure has to be used [116]. 

Since some understanding of the theory behind multigrid is necessary in order to use 
it effectively, a brief development of the Full Approximation Storage (FAS) multigrid for 
a non-linear problem is presented. Consider the problem 

L h U h = f h , (5.10) 

where L h is a non-linear operator on a grid, g h , with spacing h. The forcing function, /, 
is known and U h is the solution to the problem on the grid with spacing h. Taking u h 
as an approximation to U h with an error 

V h = U h -u h , (5.11) 

Equation (5.10) can be written as 

Wu* + V'*)=/\ (5.12) 

L h u h is subtracted from both sides of equation (5.12) to give: 

L h (u h + V //l ) - L h (u h ) = f h ~ L h (u h y (5.13) 
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If the terms are smooth, they can be represented on a coarser grid, g 2h with spacing 2h. 
The grid g* h is formed by deleting every other point in g h ; therefore, g 2h t g h . Points are 
eliminated from g 2h to form g 4h and so forth to form g 8h , g I6h , etc. Each subsequent grid 
is a subset of the previous grid, which places compatibility constraints on the number of 
grid points in each direction. Written on the coarse-grid, g 2h , equation (5.13) becomes 
L 2h (/2 V + V 2h^_ L 2 h (/2V) = /f (/ A - I V) , (5.14) 

or 

L 2h (u 2h> ) = f 2h , (5.15) 

where 

j 2 h = j 2 h ( jh _ L v) + L 2h (ll h u h ) , (5. 16) 

and l\ h is the restriction operator. 

Since equation 5.15 is on a coarser grid than equation 5.10, the numerical solution for 
u* h is much cheaper to obtain because fewer points are involved. Note that the operator 
used on the coarse-grid has the same form as the fine-grid operator, the grid spacing (h 
and 2 h) being the only difference. Once the values of u 2h are obtained, the fine-grid 
iterative solution is updated using the following equation: 

(uM =(u h ) +l!} h \u 2h -I 2 h h (u h ) 1 (5.17) 

V / New \ / Old I V / Old! 

where is the prolongation operator. 

It should be emphasized that the prolongated term on the right-hand side of equation 
(5.14) is the correction to be applied to the fine-grid solution. Examination of this term 
shows that the solution on the coarse-grid is actually a solution to the originally posed 
problem, which allows the use of the fine-grid boundary conditions on all the coarse-grids 
as well. In the developed algorithm, the non-linear FAS scheme utilizes the same operator 
on all the grid levels. This of course simplifies the programming of the multigrid scheme. 
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A grid with spacing 4 h can then be used to find corrections to the “solution” of the 
problem on the grid with spacing 2 h. Successively coarser grids may be used until a 
grid is reached which is so coarse-that a direct solution may be used (or a nearly exact 
solution with only a small number of iteration sweeps). The correction from the coarsest 
grid is then used to correct the correction on the next finer grid; and this is continued 
through successively finer grids until the finest level is reached and the approximate 
solution is updated. 

The usefulness of corrections obtained on a coarser grid is dependent on the smooth- 
ness of the fine-grid error passed to the coarse-grid. Hence, it is absolutely necessary that 
the high-frequency components of the error on the fine-grid be minimized, if not com- 
pletely eliminated. It is the responsibility of the smoother (modified Runge-Kutta and 
Predictor-Corrector Schemes) to damp the high frequency components of the error. The 
removal of the low-frequency components of the error is unimportant for all but the coars- 
est grid since these frequencies can be resolved on the coarser grids where they become 
high frequencies. If the high frequencies are not damped, then the restriction operator will 
pass aliased information to the coarser grid and the entire multigrid scheme will cease to 
converge, [52], Obviously, the choice of the smoother is critical to the proper functioning 
of multigrid. Thus the choice of the modified Runge-Kutta coefficients should be tailored 
to improve the damping properties rather than for a slight increase in the maximum CFL 
number, as discussed previously. Failure of the 1-2 Scheme to damp the high frequency 
errors, over a wide range of CFL number (Fig. 4.8), disqualified that scheme for use in 
multigrid applications. 

The cycle of work performed starting on the finest grid, successively treating the 
coarser grids, and then returning to the finest grid is called one multigrid cycle. The cycles 
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are repeated until sufficient convergence is obtained on the finest grid. In the present 
study, fixed cycles known as V- and W-cycles are used, and are given in Appendix D. 

The restriction operator has two forms. One form is used to restrict the dependent 
variables, if (tt*); i.e., the flow quantities p , pu , pv,pw, and e. For these, the volume 
weighted average of the values of the function at mid-cells of the eight fine-grid cells, 
contained in a coarse-grid cell, is used to set the value on the coarse-grid and is given 
by [85] : 

(5-18) 

*=i J2 Voi h t 

The other form of the restriction operator is for the restriction of residuals, if [L h (u h )\ . 
Following Cannizzaro et al [85], a simple summation of the residuals over the eight 
fine-grid cells composing the coarse-grid cell is performed such that; 

L 2h (u 2/l ) = if [. L h (u*)] = Y, (5.19) 

«=i 

The restriction operations are performed for all interior points. At the inflow/outflow 
boundaries, only the values of the functions are restricted, with no residual restriction. 
The residual values are frozen to the fine-grid values and are not updated on the coarse- 
grids. On wall surfaces, the same boundary conditions are used for all the grids. 

The prolongation operation used in the current work was a tri-linear interpolation, in 
the computational space, of the corrections at the eight coarse-grid cells adjacent to the 
fine-grid-mid-cell. A practical approach to the coding of a multigrid scheme in Fortran 
V is presented elsewhere [54, 117]. 

A constant coefficient implicit corrector smoother was used to remove high frequency 
errors from the coarse-grid corrections before they were applied to the fine-grid. For the 
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test cases investigated in this study, the correction smoothing procedure did not enhance 
the rate of convergence, but it should pay off for high speed flows, [114]. This operator 
is identical to that used to smooth the residual in the previous section given by eq. 5.3. 
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CHAPTER 6 
MULTI-BLOCK 

The demand to resolve the fine details of physical flows has challenged many 
researchers to find new and efficient computational tools. This challenge has pushed 
the development of new computer architectures and numerical techniques which permit 
the solving of complex, realistic geometries and configurations for compressible flow 
problems. Many problems in engineering are solved using body fitted coordinate systems. 
However, many aerodynamics designs are often quite complex (geometrically), and quite 
often, generating a single, body fitted grid for realistic three-dimensional geometries is a 
difficult task to perform; for some configurations it is almost impossible. 

In the present study, a multi-block strategy is employed to allow greater geometric 
flexibility on structured grids. The multi-block strategy (multizone) has a number of 
advantages. It alleviates the problem of grid generation for complex configurations [85]. 
Different types of governing equations can be used on different domains [75]. Multi- 
block can even allow the use of different numerical techniques and grid topologies on each 
block [54], Multiblock also requires less memory if each zone is solved independendy. 

Several grid methodologies such as, overlaid grids [83], patched grid [118, 119], and 
blocked grid [85, 120] can be applied to simplify the grid generation, as well as provide 
geometric flexibility, and mesh refinement. In the present study, blocked grids have been 
used because the flow properties are conserved automatically across the block interface. 
The result section will show that this allows discontinuities, within the computational 
domain, to move freely across the block interfaces. The multi-block strategy along with 
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the interaction between multi-block, multigrid, and time integration schemes will be 
discussed in the this chapter. The different boundary conditions applied on the block 
faces will also be discussed. 

6.1 Multi-Block Strategy 

A multi-block strategy is used to allow greater geometric flexibility. The solution 
domain is divided into multiple zones (blocks) and the grid for each zone is then generated. 
Each block within the computational domain is treated as a three-dimensional box. Each 
block can have a different grid topology. Different grid topologies are often better suited 
for a particular flow component or configuration within the computational domain. If the 
blocks, and block grid topologies are chosen appropriately, the difficulty of generating a 
boundary fitted grid can be reduced. Also, the placement and control of wall boundary 
conditions are more flexible. The trade-off is the computational overhead required 
for communication between the multiple blocks across their respective intersections 
(interfaces). In reality, the lagging of communication across the interfaces can slow 
convergence. 

Numerical treatment of grid interfaces is of paramount importance for algorithms that 
employ different grids within the computational domain. Interface boundary conditions, 
if not handled properly, can cause the numerical solution to degrade at interfaces [121], 
In the present multi-block implementation, the grid in adjacent blocks, connected across 
an interface, is assumed to have C° continuity. The grid lines at the block interface are 
continuous but the slopes are not necessarily continuous. Having C° continuity greatly 
simplifies the handling of the boundary conditions across the interface, and avoids the 
necessity of spatial interpolation of the data when loading ghost cells at the block interface. 
It also ensures the accuracy and conservation of flow properties across each interface. 
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Figure 6.1 Schematic of a Block Face with a Generic 
“Patch”, Accommodating Multiple Boundary Conditions 


This allows discontinuities within the computational domain to move freely across these 
interfaces as will be shown in the result section. 

In the early stages of the present work, the multi-block strategy [85], used a 
homogeneous boundary condition for any given face of the block. (The entire face of the 
block had to be a wall, an inflow/outflow boundary, or an interface with another block.) 
That limitation has been relaxed to allow multiple boundary conditions per face [42, 
117]. The face of each block can be divided into rectangular patches, where each patch 
can u tiliz e a different type of boundary condition, as shown in Fig. 6.1, thus increasing 
the flexibility of the code for handling complex three-dimensional configurations with 
different boundary conditions. The boundary conditions for each patch on each block 
face can be specified in an input file to the algorithm, and they are not “hard-wired” in 
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the source code. 



On the block faces that have either a wall or an inflow/outflow boundary condition, 
standard boundary conditions are used, as will be discussed later. On faces that are 
interfaces, a special interface routine presets the values in two layers of ghost cells 
(normal to the face) equal to the latest values in the coincident interior cells in the 
adjacent blocks. The updates of the interface ghost cells are performed before each 
iteration in a given block. The iteration on each block can then proceed without the need 
for further information from adjacent blocks. Hence, each block is thickened by two 
ghost cell shells which carry the solution from adjoining blocks, as boundary conditions, 
into the computationally active block. 

6.2 Multi-Block and Multigrid 

There are two possible strategies for the implementation of multigrid with a multi- 
block grid structure. Either multigrid inside of multi-block, or multi-block inside of 
multigrid can be used. The first strategy implies that a complete multigrid cycle (or 
cycles) will be performed for a given block. Subsequently, computation moves to the next 
block and so forth until all the blocks are complete. This strategy can be advantageous 
since it allows the flexibility of different numbers and/or types of multigrid cycles for 
different blocks, and they can be adjusted to speed convergence in slowly converging 
blocks (assuming only steady-state results are sought). Unfortunately, communication 
between the blocks is reduced which slows convergence. 

With the multi-block inside of multigrid strategy, the multi-block structure is just 
a way to update all the points on grid h in the multigrid cycle. Then a restriction is 
performed on all the blocks and the multigrid process is continued on each block for grid 
2h. This continues for each of the multigrid grids, and allows communication between 
the coarse-grids in the multigrid cycle, through updates of the interface conditions. It 
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or invent some interface boundary condition on the coarse-grids. This method can also 
reproduce the convergence history of a single block solution using an explicit algorithm, 
and can be used to validate the multi-block logic. The multi-block inside of multigrid 
strategy was used in the present work. A schematic of the strategy is shown in Fig. 6.2. 

6.3 Boundary Conditions 

When solving computational fluid dynamic problems, several types of boundaries can 
be encountered. These boundaries can be real boundaries or artificial boundaries. The real 
boundaries can be simple solid or porous surfaces or complex wing-body junctures, while 
artificial boundaries can be far field boundaries or symmetry planes. These boundaries 
are the link by which the computational domain senses the rest of the “universe”. They 
drive the solution in the computational domain. An inappropriate boundary condition or 
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boundary procedure can have a destabilizing effect on the numerical solution. It can be the 
difference between fast convergence to steady-state or no convergence and instabilities. 
Thus it is of paramount importance, when solving fluid dynamic problems numerically, 
to select and implement appropriate boundary conditions and boundary procedures. 

In the present study, each block within the computational domain is treated as a 
three-dimensional box. This box has six faces and it is on these faces that the bound- 
ary conditions have to be applied. Several types of boundary conditions have been 
incorporated in the developed algorithm corresponding to the different test cases investi- 
gated. These boundaries are inviscid/viscous solid walls, symmetry planes, inflow/outflow 
boundaries, and the interface between blocks. To facilitate the treatment of the boundary 
conditions and the evaluation of the fluxes at the boundaries, two layers of ghost cells 
(virtual/phantom cells) are used at the boundaries, as shown in Fig. 6.3. The ghost 
cells have their own volume and directed areas, and are similar to any other cell in the 
computational domain, except they are not updated during the block computations. The 
memory allocation for each of the blocks is increased by two planes on each of the six 
faces. A description of the various types of boundary conditions and their implementation 
in the numerical algorithm will be discussed in sections which follow. 

Solid Boundary When solving the Euler Equations, the boundary conditions to be 
applied at solid boundaries are flow tangency conditions. The velocity component normal 
to the solid boundary is set equal to zero. The pressure is extrapolated linearly from inside 
the computational domain to the wall. One only needs the pressure at the wall to compute 
the inviscid flux [122], More complex boundary conditions exist in the literature but the 
current approach has been robust and produced accurate results for the test cases studied. 

When solving the Navier-Stokes equations, we enforce the no slip and the no injection 
boundary condition at the solid boundary. The pressure is extrapolated from inside the 
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Ghost Cells 



Figure 6.3 Schematic of a Plane in the Computational Domain 
computational domain. The treatment of the inviscid part of the flux is similar to the 
treatment pursued for computing the inviscid flux in the Euler Equations. In all the 
viscous test cases investigated, it was assumed that heat transfer between the fluid and 
the solid boundary was negligible; an assumption of an adiabatic wall was made. The 
thermodynamic properties in the ghost cells, for physical boundaries, are set equal to 
the properties in the cells adjacent to the boundary. The velocity in the ghost cell is 
computed by requiring the average of the velocity in the ghost cell and the cell adjacent 
to the boundary to be equal to zero on solid boundaries. 

Symmetry Plane The third type of boundary condition encountered in this study is the 
symmetry plane boundary condition. The values in the ghost cells are set to be the mirror 
image of the interior cells at the symmetry plane. The evaluation of the flux is the same 
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as the evaluation of flux at any point in the computation domain. The only difference is 
that the symmetry plane flux utilizes information stored in the ghost cells. 


Inflow/Outflow Boundary. A non-reflective type of boundary condition at the far field 
is essential to minimize the reflection of non-physical outgoing disturbances. Thus a 
characteristic non-reflective type of boundary condition is used to compute the flow 
variables in the ghost cells at the inflow/outflow boundaries. This type of boundary 
condition is based on characteristic variables and the assumption that the flow is steady 
and locally homentropic at the boundary. The procedure for implementing the boundary 
condition was developed by Thomas and Salas [123] for two-dimensional flows. The 
derivation and application of the characteristic boundary conditions for three-dimensional 


flows are discussed in detail elsewhere [100]. 

In applying the characteristic boundary conditions for the developed algorithm, one 
can proceed by computing the Mach number normal to the boundary, at the first interior 
cell in the computational domain. That Mach number is used to determine the nature 
of the flow; subsonic or supersonic flow. If the flow is subsonic, the two Riemann 
invariants R + , and R~, are computed as: 


and 


R + = q,nt + 
R = Qref 


2 

r a t nti 

7 - 1 
2 

T a re / 1 

7 “ 1 


(6.1) 


where a and q are the speed of sound and the contravarient velocity normal to the 


boundary. The subscript int and ref indicate the first interior cell to the boundary and at 
reference conditions, respectively. Here, q re f, and q int are given by: 
q r ef = u reflx "b v ref^y "b ^ re y/ z . 


<7»nt — u tnt^x "b v t nt^y "b , 

In 


( 6 . 2 ) 


/„ = 


+ /J + /2 


where / = £,??, (, and subscript n = x, y, z 
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We can add and subtract the two Riemann invariants to compute the local velocity normal 
to the wall, q g , and the speed of sound a g such that 


q g =0.5(f? + -R~) , 
a g = 0.25(7 -1){R + -R~) 


(6.3) 


If the computed q g , normal velocity, is positive, the boundary is a subsonic outflow. 
Thus the entropy and the tangential velocity are carried outside the computational domain 
by the outgoing characteristic waves. The Cartesian velocity components and entropy in 
the ghost cells are then computed as follows 

U g ~ u int "i" (<?n — 

V„ — V{ n i + (q n Qint)^yi 

0 6 . 4 ) 

Wg = V in t + (^n — <7i nt)lzi 
s * = Pint / P' nt 

where s* is an entropy related function. 

If q gt is computed to be negative, the boundary is a subsonic inflow. Thus the entropy 
and the tangential velocity are carried inside the computational domain by the incoming 
characteristic waves. The Cartesian velocity components and entropy in the ghost cells 
are then computed as 

tig — tt-ref 9rc/)^xi 

V g — V re f + (<7n — 9re/)^yi 

( 6 . 5 ) 

Wg — Vi T ef "4" (^n 9rc/)(zi 
5 = Pref / P re f- 
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Figure 6.4 Inflow and Outflow Boundary Conditions 
The density, pressure, and the total energy in the ghost cells are computed for both the 
subsonic inflow and outflow as: 


and 



p g = 


E g = 


a ~ Pg 




( 6 . 6 ) 


For supersonic inflow, all characteristic waves are entering the computational domain 
and the flow variables in the ghost cells are set to the freestream values. In the case 
of supersonic outflow, all characteristic waves are leaving the computational domain 
and the values in the ghost cells are extrapolated from the interior. Second-order 
extrapolation was always implemented. Figure 6.4 shows the four different scenarios 
for the inflow/outflow boundary. 
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It should be emphasized that the characteristic boundary conditions are valid for 
steady, locally homentropic flow. For viscous dominated regions, a simple, robust and 
reasonable procedure is to extrapolate all the variables at the downstream boundary from 
the interior. 

Interface Between Blocks This type of boundary condition only arises in a multi-block 
domain at the interface between various blocks. This boundary is not a physical boundary 
and it is of paramount importance to treat the computational cells at the interface with 
the highest level of care, in order to ensure accurate transfer of data from one block to 
another. The interface should be transparent to the flow of information across it. 

Two ‘ghost’ cells are used to pass all the necessary information from a neighboring 
block to a cell at the face of a given block without degradation of accuracy at the interface. 
A special interface subroutine was used to load the data from the internal cells in the 
neighboring blocks into the proper ghost cells. This process was performed after the other 
boundary conditions were enforced, and between each sub-iteration on a given multigrid 
level. The flux at the interface is calculated in the same way as it is calculated at any 
other point within the computational domain. The only difference is that the interface 
plane utilizes information from its neighboring block. It should be noted that the present 
algorithm has the capability for a block to interface with itself. 
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CHAPTER 7 
TEST CASES 


7.1 Background 

The developed algorithm was utilized to compute several test cases of interest. 
The main objective of the computations were to validate, demonstrate, and assess the 
predictive capability of the algorithm. In the present chapter, results of the computation 
of corner flow, plume flow, laminar and turbulent flow over a flat plate, an ONERA M6 
wing, and the unsteady impingement of a jet on a ground plane are reported. Each test 
case was computed in order to verify a certain aspect of the developed code, as will be 
shown in the following sections. 

The first test case considered was an inviscid, three-dimensional supersonic flow 
through a comer of intersecting ramps. That flow has been investigated both experimen- 
tally and numerically by several researchers, [124-127]. The flow field encountered is 
highly complex and is dominated by complicated shock structures, shock interactions, 
shear layers, and shock boundary-layer interactions. A schematic of the flow structure is 
given elsewhere [127], and is repeated here as Fig. 7.1. 

Typical examples of such flows are supersonic/hypersonic inlets and wing body 
junctures. In the case of supersonic and hypersonic inlets, the associated shock structure 
is of paramount importance, since that flow is convected through the combustion chamber 
and can lead to non-uniform combustion. Therefore, the development of efficient 
algorithms which can predict accurately the flow field, while consuming reasonable levels 
of computer resources is highly desirable. Clearly the present numerical scheme should 
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Figure 7.1 A Schematic Diagram of a Cross Plane for Supersonic Comer Flow. 


be capable of accurately resolving the complicated shock structure on a reasonably fine 
mesh. Even in its simplest form, (inviscid formulation), the flow problem considered is of 
significant physical importance, and represented a challenge to the numerical algorithm. 


The second test case computed was the non-axisymmetric jet exhaust plume. Several 
studies have been conducted, [75, 128, 129], to investigate the performance of non- 
axisymmetric jet exhaust plumes. The motivation for the non-axisymmetric jet flow 
studies has been the fact that they can provide comparable or superior level flight 
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performance to axisymmetric nozzles, [129]. The rectangular exit nozzle provides the jet 
engine with improved thrust vectoring and thrust reversing capability, which can enable 
a substantial reduction in landing and take off distance requirements while enhancing 
the maneuverability in flight. An additional advantage of the rectangular exit over the 
axisymmetric nozzle is that the design and structure of the rectangular nozzle is simpler 
and lighter for the same capabilities. 

Propulsive nozzles of this type have been considered for a wide range of possible 
applications in airplanes such as the Advance Tactical Fighter (ATF), the short takeoff 
and landing (STOL) Eagle F-15 (F-15B) and the National Aerospace Plane (NASP). 

A better understanding of the type of flow surrounding the nozzle and plume 
region is needed because of its influence on the sonic boom signature and aerodynamic 
performance, [128]. 

The non-axisymmetric jet exhaust plume flow is complex physically and represents a 
challenging test case for the algorithm developed in this study. Expected exhaust plume 
mixing (with the surrounding flow) will occur at transonic conditions, accompanied by 
large embedded supersonic flow regions and possibly complicated shock structures, shock 
interactions, slip lines and shear layers. A schematic of the flow field showing the shock 
wave, expansion fan, and slip line are demonstrated in Fig. 7.2. 

Laminar and turbulent flows over a flat plate were computed to verify the implemen- 
tation of the viscous terms and turbulence model (Baldwin-Lomax algebraic-turbulence 
model). The turbulent flow over an ONER A M6 wing was also computed in order 
to compare the performance of the developed algorithm with other three-dimensional 
state-of-the-art computer codes. 

The ONERA M6 wing is a basic three-dimensional configuration. The flow field 
around the wing contains a wealth of aerodynamic features such as; shock waves. 
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spanwise variations in boundary layers, and the interaction between shock waves and 
the boundary layer. The wing has been developed for the experimental support of 
three-dimensional, transonic and subsonic flow fields and has been used extensively 
as a benchmark case to gauge the accuracy and performance of newly developed 
computational codes. Both Euler and Navier-Stokes solutions have been reported by 
several researchers [130-133], and there is also a large data base of experimental 
pressure distributions available, [134], Hence the ONERA M6 wing has been selected to 
validate the capability of the developed algorithm to compute a truly three-dimensional 
turbulent-flow, and to access the performance of the developed algorithm with other 
three-dimensional, state-of-the-art computer codes. 

The flow of a jet impinging on a ground plane was the last test case conducted to 
validate the developed algorithm. Impinging jets, with and without cross flow, occur 
in a wide variety of engineering applications (cooling, heating, and drying of a variety 
of industrial products, tempering of glass, cooling of turbine blades, paint spraying, and 
welding, are just a few). Particular attention has been devoted to the problem, because 
of its application to V/STOL (vertical/short takeoff and landing) aircraft [135-147], A 
comprehensive review, of previous experimental and numerical work, for an isolated jet 
impinging on the ground plane is given by Jalamam [137]. An up-to-date review of 
impinging jets in cross flow is given by Bray and Knowles [144], The associated flow 
field of impinging jets, is highly complex, and unsteady. The flow field contains highly 
sheared layers, vortical regions, an impingement zone, free jet, and a wall jet. The jet 
impingement problem is a demanding test case. The analysis of such a complicated 
flow field, requires the solution of the three-dimensional, time-dependent, compressible, 
Navier-Stokes equations, to be capable of accurately resolving the entire flow field. 
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A Schematic of the flow field of an isolated jet impinging on a ground plane is shown 
in Fig. 7.3. The numerical simulation of such a complicated flow field will demonstrate 
the predictive capability of the developed algorithm, to resolve complex unsteady flows. 




Figure 7.3 Schematic of an Isolated Jet Impinging on a Ground Plane. 


When a jet impinges on a ground plane a wall jet is formed. The wall jet flows 
radially outwards from the zone of impingement. In the presence of a crossflow, the 
wall jet flowing radially outward is opposed by the crossflow (free stream) and rolls up 
into a horseshoe-ground vortex as shown in Fig. 7.4. The ground vortex is the primary 
mechanism for hot gas ingestion, creating dust clouds, and causing lift loss for the 
V/STOL aircraft. 

There is a significant amount of scatter in the database available for jet impinging 
studies [144], Extra numerical and experimental research work is highly recommended 
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to complete the development cycle of V/STOL aircraft. Developing a tool that can be 
used to simulate numerically the flow field of a jet impinging on a ground plane, with 
and without cross flow is highly desirable. 



Figure 7.4 Schematic of a Jet Impinging on a Ground Plane in Presence of Crossflow 

Even in its simplest form, the flow problems to be considered in the proposed work 
are of significant physical complexity. The computed flow fields will contain most of the 
rich features of fluid mechanics (shocks, rarefaction waves, shear layers, etc.). Clearly 
the numerical scheme should be capable of resolving all of these features accurately on 
a reasonably fine mesh. The developed algorithm must prove to be robust and reliable 
to provide the user with the necessary confidence in making design decisions based on 
the results obtained from this scheme. 
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7.2 Corner Flow 


The flow through a rectangular channel configuration with two compression ramps 
forming a compression comer about which the channel is symmetric was the first test 
case computed. A schematic of the compression comer is shown in Fig. 7.5. The back 
and bottom walls of the channel have converging ramps, each with inclinations of 9.5°. 
A supersonic inlet flow Mach number of 3.0 was used for the first test case. The 
2-2 Predictor-Corrector, Explicit Scheme was employed to compute the flow field. Mach 
contours on the two side- walls, and the exit plane are shown in Fig. 7.6, while the 
pressure contours are shown in Fig. 7.7. The Mach and pressure contours show clearly 
that three shock surfaces are generated. Two of the shock surfaces are two-dimensional 



94 



wedge-flow shocks, which can be verified using a two-dimensional analysis based on 
the Mach number normal to the leading edge of the wedge. The third surface is formed 
where the two wedge shocks coalesce to form a three-dimensional flow region which 
is shaped like a cone. This can be seen in Fig. 7.6, where the Mach line contours are 
shown on the back and bottom walls, and on the exit plane of the channel. The positions 
of the wedge shocks are apparent on the back and bottom walls as regions of highly 
concentrated Mach lines perpendicular to the flow direction. Also, on these two walls 
the edges of the cone shaped surface can be seen. On the exit plane, four flow regions are 
present. In the upper right corner is free stream, which is one-dimensional flow. From 
the middle of the plane to the lower left comer (where the flow is three-dimensional), 
the bottom of the cone surface appears as a partial disc. The two wedge shock planes 
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can be seen in the upper left corner and the lower right comer of the exit plane. Note 
that since the geometry of the channel is symmetric about the compression comer (the 
one joining the back and bottom walls to the exit plane of the channel), the steady-state 
flow field is symmetric about the corner. 

A schematic of the different flow zones, at the exit plane is discussed by Marconi 
[127], and is shown in Fig. 7.1. The wall shock, comer shock, two-dimensional flow 
regions, can be identified, along with the one-dimensional free stream flow, and the 
three-dimensional regions in Fig. 7.1. A qualitative comparison between the computed 
results. Fig. 7.6, and Fig. 7.7, and the schematic diagram of the flow field [127], clearly 
shows that the present algorithm captured all four zones in the flow field accurately, and 
resolved the complicated shock structures. The triple points, where the one-dimensional 
free stream region meets, the two-dimensional wedge flow and the three-dimensional flow 
region, have been captured by the present computation, as shown in Fig. 7.6 and 7.7. 

A comparison between the computed results, the two-dimensional results of Marconi 
[127], and Kutler [126], and the experimental data of Charwat and Redekeopp [124] is 
shown in Fig. 7.8. Marconi [127] used a shock-fitting code to obtain his results, while 
Kutler [126] used a second-order central-difference shock capturing scheme. The present 
results are in good agreement with the numerical and experimental results. The present 
numerical results, as well as Kutler’s [126] results, suffer from over- and under-shoots 
at the shock wave, as expected. It should be emphasized that these oscillations are due 
to the fact that the shock wave is cutting obliquely across the grid lines. Pressure 
oscillations near the shock wave are evident in Kutler’s results. The reference pressure, 
Pb, in Fig. 7.8, is the pressure value on the side wall in the two-dimensional flow region. 
Y 0 is the point where the two compression ramps intersect, corresponding to the location 
of Pb, as shown in Fig. 7.5. 
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To compare the performance of the different time integration schemes, the inviscid 
corner flow was recomputed on a 33x33x33 grid. Comparison between the convergence 
history to achieve steady-state for the explicit two-stage predictor-corrector schemes 
(1-2 and 2-2 Schemes) and the four-stage « = 0 explicit scheme are shown in Fig. 7.9. 
It should be emphasized that the predictor-corrector schemes need only about half the 
CPU time per iteration, compared to the four-stage, explicit time-stepping schemes. 
The convergence history to steady-state for the four-stage explicit scheme is shown in 
Fig. 7.10. By comparing Fig. 7.9 and Fig. 7.10, it can be seen that the best performance 



Surface Pressure Distributions of Comer Flow. 


98 





was achieved with the 1-2 Scheme for this particular problem. The main drawback of 
the 1-2 Scheme is its incompatibility with multigrid acceleration techniques as shown 
in Fig. 7.11. This is in agreement with the stability analysis study reported earlier in 
chapter four. The second drawback was that the convergence rate of the 1-2 Scheme 
deteriorates on fine (viscous) grids. 

To validate the multi-block capability of the algorithm, and to ensure that the block 
interface is transparent to the numerical scheme, the inviscid corner flow was recomputed 
using eight blocks. The previous single block grid was used as a starting point to generate 
the grids in the eight blocks by dividing the domain in half in each of the three coordinate 
directions (Fig. 7.12). The Mach contours on the two side walls, and the exit plane for the 
eight block calculation are shown in Fig. 7.13, while the pressure contours are shown in 
Fig. 7.14. Results of the multi-block calculation reproduce the results obtained with the 
single block calculation identically. As shown in Figs. 7.13 and 7.14, the block interface 
is transparent to the developed algorithm. 

The convergence histories for the single-block and the eight-block calculations are 
shown in Fig. 7.15. Notice that the convergence rate for the two-step explicit scheme 
shows little degradation for the multiple block calculation. This is due primarily to 
the choice of the multi-block-inside-of-multigrid strategy which allows communication 
between the coarse-grids in the multigrid scheme. The small differences between the 
two convergence histories are due to the differences in overhead due to multi-block data 
transfer. 
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Figure 7.9 Comparison of Convergence Histories of 1-2, 2-2, Four-Stage, 

K = 0 Schemes Calculations of 9.5° Compression Corner Flow (M in iet = 3.0) 
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Figure 7.10 Comparison of Convergence Histories of the Four-Stage 
Schemes, Calculations of 9.5° Compression Comer Flow = 3.0) 
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Figure 7.11 Comparison of Convergence Histories of 1-2 Scheme With and Without 
Multigrid Acceleration, Calculations of 9.5° Compression Comer Flow (M inlet = 3.0) 
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Figure 7.12 Schematic of Grid for Eight Block 
Calculation for Flow through a 9.5° Compression Comer 
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Figure 7.15 Comparison of Convergence Histories of Single Block and Eight 
Block Calculations of 9.5° Compression Comer Flow ( M, n i et = 3.0) 
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7.3 Non-Axisymmetric Jet Exhaust Plume 

Non-axisymmetric jet exhaust plume flows generally contain complex, multiple 
shocks and strong contact discontinuities or slip lines. The flow complexity is driven 
by both geometric complexities and Mach number, but is also influenced by static 
temperature, and static pressure ratios between the jet and the free-stream flow. The 
non-axisymmetric jet exhaust plume flow represents a challenging test case to the 
developed algorithm. The exhaust plume mixing with the surrounding flow will occur 
at transonic conditions, accompanied by large embedded supersonic flow regions and 
possibly complicated shock structures, shock interactions, slip lines and shear layers. 

The main emphasis of this test case, besides understanding the flow physics, was 
to assess the performance of the Roe flux-differencing scheme and van Leer flux-vector 
splitting scheme, for the prediction of jet exhaust plume flows. A test case was also 
conducted to determine which type of extrapolation (conservative or primitive) should 
be ut iliz ed to evaluate the cell interface flux. This test case was used to verify the 
non-homogeneous boundary condition capability of the code as well. 

To isolate the problems that Roe’s flux-difference splitting scheme and van Leer’s 
flux-vector splitting scheme have with shocks and contact discontinuities from problems 
with geometric complexities, a simple, pseudo-two-dimensional test case was considered. 
The test case was a flow from an infinitely wide nozzle of height 1.0. The jet Mach 
number was taken to be 1.5 and the free-stream Mach number was 2.5. The ratio of the 
jet static pressure to the free-stream static pressure (Pjet/Po o) was 3.5 and the ratio of 
the jet static temperature to the free-stream static temperature ( T^tlT^ ) was 3.0. 

Although three-dimensional calculations were performed using both Roe’s flux- 
difference splitting and van Leer’s flux-vector splitting schemes, they were compared 
with a two-dimensional, shock-fitting method [148]. The flow field is symmetric about 
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the x-plane and thus only half the flow field is computed numerically. Figure 7.16 is a 
schematic of the upper half of the flow field showing the shock wave, expansion fan, 
and slip line. A partial view of the grid used is shown in Fig. 7.17. When generating 
the grid, special attention was devoted to the alignment of the grid lines, as much as 
possible, with the shock wave and slip line. 

y 


■* 



Figure 7.16 Schematic of the Computed Flow Field 
for a Pseudo-Two-Dimensional Jet Exhaust Plume 


107 





io.o r 



Figure 7.17 Partial View of the Two-Dimensional Jet Exhaust Plume Grid. 


The Mach and pressure line contours for the Roe’s flux-diferencing splitting scheme 
are shown in Fig. 7.18, while the Mach and pressure line contours for the van Leer flux- 
vector splitting scheme are shown in Fig. 7.19. Both schemes employed the 1-2 explicit 
scheme. No flux limiters were used for the results displayed in Fig. 7.18 and Fig. 7.19. 

The results obtained by the Roe flux-difference splitting and van Leer flux-vector 
splitting schemes are comparable. The main difference is that the van Leer flux-vector 
splitting scheme smears the slip line slightly. This behavior is expected since van Leers s 
scheme ignores the linear waves (entropy and shear layer), and does not have a mechanism 
to resolve the slip lines accurately [107], while Roe’s flux-difference splitting is based 
on the idea of solving a Riemann problem at each cell interface. 
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a. Mach Contours. 
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a. Pressure Contours. 

Figure 7.18 Non-Axisymmetric Jet Exhaust Plume Calculations Utilizing 
Roe’s Scheme, Moo = 2.5, M je t — 1.5, Pjet = 3.5, and Xj e j/Too — 3.0. 
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a. Mach Contours. 
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b. Pressure Contours. 

Figure 7.19 Non-Axisymmetric Jet Exhaust Plume Calculations Utilizing 
van Leer’s Scheme, A/oc = 2.5, M ]e t = 1.5, Pjet/Poc ~ 3.5 , and Tj e t/T 0 0 — 3.0. 
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To allow de tail ed comparisons between the Roe, and van Leer schemes and the 
shock-fitting method of Salas [148], Mach number versus dimensionless height, y/h, 
was examined at the various downstream (x/h) locations. These plots are presented 
in Figs. 7.20a-f. The vertical grid spacing for each of the calculations is represented 
symbolically on the right of the plot. The grid spacing for the shock-fitting algorithm is 
shown as the left column. The shock-fitting algorithm generates its own grid, depending 
on the test case. The van Leer, flux-vector splitting scheme and the Roe flux-difference 
splitting scheme were computed using the same grid which was shown in Fig. 7.17. The 
spacing in the y-direction for both schemes is shown along the legend for consistency. 

In Fig. 7.20a, the Mach number distribution at x/h ~1.0 is presented. At this 
location, the expansion fan is just reaching the lower wall. Along the expansion fan, 
a region of uniform flow is present, and a contact discontinuity is observed at y / h «1.2. 
Notice that both the van Leer and Roe schemes smear this feature equally. This smearing 
is a result of the slip line cutting diagonally across the grid. After another region of 
uniform flow, shock occurs at y/h ~ 1 .65. Smearing and overshoots occur with both the 
van Leer and Roe schemes due to the skewed grid relative to that shock wave. 

At x/h «2.0 (Fig. 7.20b), the expansion fan has reflected off the bottom wall and 
is moving out toward the slip line. Again, both the slip line and shock are smeared and 
predicted with overshoots by the two upwind schemes. At x/h «2.5 (Fig. 7.20c), the 
expansion fan has moved further out toward the slip line. At x/h ~3.0 (Fig. 7.20d), the 
expansion fan is split by the slip line, as is apparent from the shock-fitting solution. This 
is not obvious in either the van Leer or Roe schemes, due to smearing and overshoots. 
At x/h «5.0 (Fig. 7.20e), the slip line clearly splits the expansion fan as indicated by all 
three methods. Notice that the shock-fitting method predicts spurious oscillations near 
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Figure 7.20 Comparison Between Roe’s Flux-Difference Splitting Scheme, and 
van Leer’s Flux- Vector Splitting Scheme and Salas’ Shock-Fitting Method 
for a Pseudo-Two-Dimensional Exhaust Plume. (Continued . . . ) 
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Figure 7.20 Comparison Between Roe’s Flux-Difference Splitting 


Scheme, and van Leer’s Flux- Vector Splitting Scheme and Salas’ 
Shock-Fitting Method for a Pseudo-Two-Dimensional Exhaust Plume. 
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the lower wall. These wiggles are more pronounced at x/h- 10 (Fig. 7.20f). At this 
station, the shock has moved out of the computational domain. 

Although overshoots and smearing are present, several important points should be 
noted from Fig. 7.20. Both the van Leer and Roe schemes do a reasonable job of 
predicting the flow and are not affected by the overshoots and smearing adversely. It 
was anticipated that the Roe scheme would do a better job of predicting the contact 
discontinuity than the van Leer scheme. This is not apparent in the present results and 
should be attributed to the ’poor’ predictions by both methods, due to the misalignment of 
the grid. Obviously, a better grid can be generated but the current grid is more indicative 
of the type of grid that would be used for more complex geometries and flow physics. 
For these complex situations, proper alignment of the grid would be impossible without 

an adaptive grid scheme. 

Evaluation of the cell interface flux requires the extrapolation of state variables to 
the cell interface. The developed algorithm offers the capability of extrapolating either 
the conservative or primitive variables. To examine the effect of utilizing either of 
the extrapolation methodologies, the four-stage k = -1, Roe’s flux-difference splitting 
scheme, with primitive and conservative extrapolation, was employed to recompute the 
two-dimensional plume flow. Local time-stepping, implicit residual smoothing and the 
multigrid acceleration technique were implemented to accelerate the rate of convergence 
to steady-state. Comparison between the computed Mach number distribution with the 
shock-fitting method of Salas [148], at various x//i-locations have been made. These 
plots are presented in Figs. 7.21a-f. The results show clearly that the primitive variable 
extrapolation renders a smoother solution across slip lines and shocks, for all x/h- 
locations. It should be mentioned that the results computed by the predictor-corrector 
explicit scheme and the four-stage k = -1 explicit scheme were identical. 


115 





(a) 


4 r 


y/h 2 


Conservative. 
Primative. 
Shock Fitting 




X 


X 


x/h ** 2.5 


_L 


_L_ 


1.7 1.9 2.1 2.3 2.5 

Mach 


(b) 

Figure 7.21 Comparisons Between Primitive and Conservative 
Extrapolations of the Roe’s Scheme with the Shock-Fitting Code of 
for a Pseudo-Two-Dimensional Exhaust Plume. (Continued . . 
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Figure 7.21 Comparisons Between Primitive and Conservative 
Extrapolations of the Roe’s Scheme with the Shock-Fitting Code of Salas 
for a Pseudo-Two-Dimensional Exhaust Plume. (Continued . . . ) 
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Figure 7.21 Comparisons Between Primitive and Conservative 
Extrapolations of the Roe’s Scheme with the Shock-Fitting 
Code of Salas for a Pseudo-Two-Dimensional Exhaust Plume. 
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To examine the effect of using a higher order method, the Roe scheme was recom- 
puted as first, second (« = -1, and « = 0), and third-order accurate (k = 1/3). Comparison 
between the computed results and the shock-fitting method of Salas [18], for the Mach 
number versus y/h at xjh « 1.0 is presented in Fig. 7.22. The first-order Roe scheme is 
highly dissipative and gives essentially useless results. The second and third-order results 
are more accurate and are nearly identical. The second-order (k = -1) generates a larger 
undershoot at the slip line, but it produces the least amount of overshoot at the shock. 
The third-order accurate (k = 1/3) method generated oscillations and produced the largest 
overshoot at the shock. For detailed comparison of the different order of extrapolations 
at all other locations, the interested reader is referred to reference [54], Comparison of 
the convergence history for all four types of extrapolations are shown in Fig. 7.23 



Figure 7.22 Comparisons Between Different Extrapolations of Roe’s Scheme 


and a Shock-Fitting Code for a Pseudo-Two-Dimensional Exhaust Plume. 
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Figure 7.23 Comparisons of Convergence History 
for Pseudo-Two-Dimensional Exhaust Plume. 

To examine the flexibility of the multi-block structure, and the non-homogeneous 
boundary conditions, the flow was recomputed on a two-block structured grid. The 
computation domain was divided into two blocks! one block had a homogenous inflow 
boundary condition, which is the jet exhaust, Mj et = 1.5), and a second block which had 
a homogenous free stream inflow, M^c = 2.0, (refer to Fig. 7.16). The two block results 
were identical to the single block results. No modifications to the code were required to 
go from the single-block plume calculation, with non-homogeneous boundary conditions 
at the inflow, to the two-block calculation. Only the input to the program was changed 
to accommodate both runs. 
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7.4 Flat Plate 


The main goal of computing the laminar, and turbulent flow over a flat plate was 
to verify the viscous capability of the developed algorithm, and to check the correctness 
of the Baldwin-Lomax algebraic-turbulence model implementation. For the laminar case 
the free stream Mach number (Moc) was 0.5, and the free stream Reynolds number 
(Reoo) was 1000 per unit length of the flat plate. The grid used for this test case was 
65x65x5 (streamwise, normal, spanwise). The normalized minimum spacing in the 
normal direction to the wall was. Ay, = 1x10^. The grid was fine enough to produce 
at least ten vertical grid points through the boundary-layer thicknesses at all locations on 
the plate. Although the flat plate problem was a two-dimensional problem, five spanwise 
planes were employed to test the multigrid acceleration technique, and to check the 
viscous terms in all three directions. The computed laminar velocity profile, and variation 
of the local skin friction coefficient along the plate are shown in Fig. 7.24. Comparison 
between the computed results and the classical Blasius boundary-layer solution , [93], 
shows excellent agreement. 

A grid refinement study was performed to determine the minimum number of points 
required to accurately resolve the laminar boundary layer. Three grids were employed 
in the grid refinement study; the first grid was a 33x49x5 with a Ay, =0.015, the 
second grid was 33x65x5 with a Ay, = 0.0075, and the third grid was 33x81x5 with 
a Ay, = 0.00375. The grid refinement was carried out by dividing the grid spacing 
in the normal direction only. Comparisons between the computed velocity profiles, and 
Blasius boundary-layer solution are shown in Fig. 7.25. The developed algorithm has 
resolved the laminar boundary layer successfully, even on the coarsest grid, where there 
were only five points in the boundary layer. Fig. 7.26 shows a comparison between the 
computed skin friction and Blasius boundary-layer solution, for all three grids. 
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Figure 7.25 Velocity Distribution Over a Flat Plate. 





The flow over the flat plate was recomputed to evaluate the effect of the implicit 
residual smoothing, and the multigrid acceleration techniques on the rate of convergence 
to steady-state. The recomputed case employed the 33x65x5 grid where, Ayi = 0.0075. 
Four sets of computations were performed. In all four sets the four-stage k = 0 explicit 
time-stepping scheme was used to drive the solution to steady-state. In the first set, 
the basic explicit algorithm was augmented with local time-stepping. In the second set 
of computations, the basic explicit scheme was coupled with local time-stepping and 
implicit, residual-smoothing. In the third set of calculations the basic explicit scheme 
used local time-stepping and the full multigrid acceleration technique. In the last set of 
calculations, the performance of the basic explicit scheme was boosted by utilizing local 
time-stepping, implicit residual smoothing and the full multigrid acceleration technique. 
Comparison of the convergence history for all four test cases are shown in Fig. 7.27. 
As shown in the figure the acceleration techniques do enhance the rate of convergence 
to steady-state. The figure demonstrates clearly that the implicit-residual smoothing is 
beneficial in reducing the total number of work units required for convergence. 

It should be noted, that the high frequency oscillations present on the fine mesh are 
slowly damped for highly skewed cells. However it was found that the convergence of 
the global flow field is improved significantly by using multigrid techniques, especially 
for complex three-dimensional grids. 

For turbulent flow over a flat plate, the free stream Mach number (Moo) was kept at 
0.5, and the free stream Reynolds number (Reoo) was set equal to l.Ox 10^ per unit length 
of the flat plate. The grid used in this case was 65x65x5, with a Ayi = 1.0 x 10 -5 . 
Computations were performed using the four-stage, k = 0, (second-order, upwind-biased) 
explicit, time-stepping scheme. Based on the previous experience gained when solving the 
laminar flat plate, local time-stepping, implicit-residual smoothing, and the full multigrid 
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acceleration technique were employed to enhance the rate of convergence to steady- 
state. Results for the turbulent flow over a flat plate are presented in Fig. 7.28. The 
computed results for the coefficient of friction agree well with the empirical formula 
C f = 0.0592 Re~ 2 , [93], especially once fully developed turbulent flow is reached. 
Comparison between the computed velocity profile and the 1/7 power law, shows fairly 
good agreement. The comparison between the computed velocity profile with the law of 
the wall also shows good agreement. 



Figure 7.27 Convergence History For Flow Over a Flat Plate. 
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Figure 7.28 Turbulent Flow Over a Flat Plate, 
(Moo = 0.5, Reoo = l.OxlO 6 per Unit Length. 


7.5 ONERA M6 Wing 

The ONERA M6 wing was selected to validate the capability of the developed algo- 
rithm to compute a truly three-dimensional turbulent-flow, and to assess the performance 
of the algorithm with other three-dimensional, state-of-the-art computer codes. 

In the present study, a 193x49x33 (streamwise, normal, spanwise) grid was em- 
ployed. The grid employed a C-0 mesh topology; C in the streamwise direction and O 
in the spanwise direction. The grid was generated by Bruce Wedan [149], and is shown 
in Fig. 7.29. Three test cases were investigated; one was a subcritical flow and two 
cases were supercritical flows. On all three test cases, a four-stage k 1 explicit time- 
stepping scheme is used to compute the flow field around the ONERA M6 wing. Local 
time-stepping, implicit residual-smoothing, and full multigrid procedures were added to 
the explicit time-stepping scheme to accelerate convergence to steady-state. Roe’s flux- 
differencing scheme was employed to compute the inviscid flux, and a second-order, 
central-difference operator was used to evaluate the viscous and diffusion terms. The 
second-order upwind scheme was selected for these computations, to avoid the use of a 
limiter when computing the supercritical test cases. Limiters halt convergence and lead 
to limit cycle solution oscillation with no apparent convergence. 

The first test case was a subcritical case, where the free stream Mach number (Moo) 
was 0.699, the angle of attack (a) is 3.06°, and the chord-based Reynolds number (Reoc) 
was 11.7 x 10 6 . The Reynolds number was based on free stream flow conditions and 
the mean aerodynamic chord. Pressure distributions along several spanwise locations 
of the ONERA M6 wing are shown in Figs. 7.30a-f. The computed results show good 
agreement with the experimental data, [134], as shown by the C p plots in Figs. 7.30a-f. 
In Fig. 7.30, r) represent the locations where the pressure distribution is displayed and 
is equal to where Y is the spanwise distance and is normalized with respect to the 
semispan B/2. 
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Figure 7.29 Partial View of C-0 Grid Topology fa - ONER A M6 Wing. 
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Figure 7.30 Pressure Distribution for ONERA M6 Wing 
Moo = 0.669, a — 3.06°, and Reoo = ll-7x 10^ 
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The second test case was a supercritical case, with Moo = 0.84, a - 3.06°, and 
Reoo = 11-7 x 10 6 . This case is considered as one of the standard test cases available 
in the literature for validation of three-dimensional CFD algorithms. The same C-0 
grid reported earlier is used for the supercritical test cases. The flow is attached over the 
whole wing. Pressure contours on the upper wing surface are shown in Fig. 7.31, where a 
Lambda shock is clearly defined. Comparison between the computed results, experimental 
data [134], and the computational results from two other numerical algorithms, CFL3D 
and TLNS3D are displayed in Fig. 7.32. CFL3D, [73], and TLNS3D, [28], are 
considered to be state-of-the-art upwind and central difference schemes, respectively. 
The results reported here were performed using the thin-layer, Reynolds-Averaged Navier- 
Stokes equation and the Baldwin-Lomax, algebraic-turbulence model. 

CFL3D is an upwind code which employs an implicit approximate-factorization time- 
marching algorithm. Roe’s flux-differencing scheme is used to compute the in viscid 
flux with k = 1/3. The viscous terms are evaluated through a second- order-accurate, 
differencing operator. TLNS3D utilizes a second-order central-difference operator for 
the spatial derivatives and an explicit, five-stage Runge-Kutta time marching scheme. 
Artificial dissipation was added to the central-difference algorithm to maintain numerical 
stability, and suppress oscillations in the vicinity of shock waves and stagnation points. 
For more detailed information about CFL3D and TLNS3D, the interested reader should 
consult the papers by the original developers of the two codes [73, 28], Figure 7.32 
shows good agreement between the present results and the numerical results obtained by 
CFL3D and TLNS3D. The computed results also agree well with the experimental data. 
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The third test case is a supercritical case, with Moo - 0.84, a - 6.06°, and 
x io 6 . Again the same gnd was used. This test case is more demanding 
than the two previous test cases, because stronger shocks develop on the upper surface 
of the wing, resulting in significant flow separation. The pressure contours on the upper 
surface are shown in Fig. 7.33. The streamlines on the upper surface of the ONER A M6 
wing are shown in Fig. 7.34, where the separation region is clearly delineated. Compari- 
son between pressure coefficients for the present computed results, previous experimental 
data [134], CFL3D [73], and TLNS3D [28], are displayed in Fig. 7.35. Large discrep- 
ancies between the computed results and the experiment exist, as shown in Fig. 7.35. As 
expected, the Baldwin-Lomax, algebraic-turbulence model cannot handle separated flow, 
and was incapable of capturing the global features of that flow. The size of the separated 
region was under-predicted, and as a result, the shock is located farther downstream. 



Figure 7.33 Pressure Contours for ONERA M6 Wing 
Moo = 0.84, and a = 6.06°, and Reoo = ll-7x 10 6 . 
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It should be emphasized that the performance of the present algorithm was comparable 
to the other two codes and the discrepancies between the computed results and experiment 
were not due to a lack of resolution of the grid, but due to limitations on the Baldwin- 
Lomax, algebraic-turbulence model. Abid, Vatsa, Johnson, and Wedan [150] reported that 
the agreement between TLNS3D results, obtained with the same grid, and experimental 
data were improved significantly when the Johnson-King turbulence model was employed. 



Figure 7.34 Wall Streamlines for ONERA M6 Wing 
Moo = 0.84, a = 6.06°, and Reoo = 11.7xl0 6 . 
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7.6 Jet In Ground Effect 


The last test case to validate the developed algorithm was a simulation of a jet imping- 
ing on a flat plane. The analysis of such a complicated flow field, requires the solution of 
the three-dimensional, time-dependent, compressible, Navier-Stokes equations, in order 
to resolve the entire flow field accurately. 

Two test cases were conducted: a single isolated-jet impinging on a ground-plane, 
and a jet impinging on a ground plane in the presence of an ambient crossflow. The jet 
impinging problem is a demanding test case. It should be emphasized that the numerical 
calculations which were performed here were devoted primarily to the development and 
verification of the algorithm for this three-dimensional, time-dependent, compressible 
Navier-Stokes case. Also, it was desired to demonstrate and assess the capability of the 
algorithm to capture the large scale phenomena associated with these unsteady flows. 

The Baldwin-Lomax, algebraic-turbulence model incorporated to date in this algo- 
rithm is not capable of resolving separated flows accurately. As demonstrated in the 
previous test case (ONERA M6 wing at 6°angle of attack), separating and non-parallel 
flow features are not modeled properly. The use of classical Reynolds-averaged turbu- 
lence models to account for such a flow field is questionable and can even mask the fluid 
mechanics to be studied [137, 151]. Developing a universal turbulence model capable 
of resolving such a complicated flow field is beyond the scope of the present work, but 
should be pursued in the future. 

Jet Impinging on a Ground Plane The numerical simulation of the flow field gen- 
erated by an isolated jet impinging on a ground plane, in ambient air, was performed 
for a jet Mach number of 0.5 and a jet Reynolds number of 19000, based on the jet 
diameter. The jet exit was assumed to be 4 jet diameters away from the ground plane, 
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and the ground-plane was assumed to extend radially outward 15 jet diameters. The flow 
properties were normalized with respect to the jet inlet conditions. 

A partial view of the grid used in the numerical computation is shown in Fig. 7.36. 
Grid lines are clustered near the solid boundary, and along the edge of the free- and 
wall-jet. The grid dimensions were 129 (normal) x 97 (radial) x 5 (circumferential) 
where the flow was assumed to be symmetric. 

A schematic of the flowfield for the isolated jet impinging on a ground plane was 
shown previously in Fig. 7.3. The flow field is characterized by distinct regions, i.e., 
free jet, shear layer, stagnation, and wall jet regions. The numerical simulation of such a 
complicated flow field can demonstrate the predictive capability of the present algorithm 
for resolving complex flows. 
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Computing a Jet Impinging on a Ground Plane. 
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The numerical simulation started with the jet entering the computational domain and 
interacting with quiescent air. Fig. 7.37a . The solution was advanced in time using a non- 
dimensional, global time step. The time step was based on the stability criterion for the 
scheme. Accelerating techniques, local time-stepping, residual smoothing and multigrid 
were not used in this test case because the flow was unsteady. Snap shots, at different time 
steps, of the Mach number and the stagnation pressure contours are shown in Fig. 7.37. 
A starting vortex developed as the free jet propagated into the computational domain, 
as shown in Figs. 7.37b-c. The free-jet impinged on the ground-plane. Fig. 7.37d, and 
is deflected to form a wall jet, as shown in Fig. 7.37e. The stagnation region creates a 
favorable pressure gradient which causes the wall jet to accelerate rapidly as it departs 
radially from the stagnation region, Figs. 7.37f-g. The primary vortex, located on top of 
the wall jet, produced a local unsteady adverse pressure gradient which caused the flow 
to separate. The wall jet, with a primary vortex located above, moves with the separation 
point radially outward as shown in Figs. 7.37f-h. From Fig. 7.37, the different flow fields 
that identify this type of flow, namely the free jet, the stagnation region, the wall jet and 
the moving separation, are predicted accurately and appear to be accurate depictions of 
the flow physics. 
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J. 






(d) 


Stagnation Pressure Contours Mach Number Contours 

Figure 7.37 Jet Impingement on a Ground Plane. M jet = 0.5, 
Pjet/Poo = 1, 7 ^- = 1, 77 = 4. Re = 19000. (Continued . . . ) 
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Figure 7.37 Jet Impingement on a Ground Plane. 


Mjet = 0.5, Pjet/Poo = 1. 


§ = 4, Re = 19000. 
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Profiles of the time-averaged velocity (tangential component), versus the normal 
wall distance off the plane are shown in Fig. 7.38. It can be observed that the favorable 
pressure gradient created by the jet impinging on the ground plane causes the wall jet 
to accelerate as it departs from the stagnation region. Fig. 7.39. As shown in Fig. 7.38, 
the time averaged velocity increased, from the stagnation region, to a maximum value 
downstream (x/r = 2.81). Further downstream the maximum value decreased and the 
boundary layer increased. The flow in the impact region (stagnation region) is dominated 
by the high pressure gradient, and the flow behaves in an almost inviscid manner. 

Unfortunately, there are no experimental or numerical data sets available to permit 
comparison of the computed time-averaged profiles. However, a flow visualization study 
conducted by Didden and Ho [138], reported the unsteady separation in the boundary layer 
of an impinging jet on a flat plane. Didden and Ho, suggested that unsteady separation was 
induced by the primary vortex and it moved downstream in the radial direction. Jalamam 
[137], computed numerically the two-dimensional time-dependent, impulsively started jet 
issuing from a plate in proximity to the ground. His results indicated the formation of a 
primary vortex, and separation of the wall jet as it moved radially outward. 

Comparison between the computed ground pressure distributions, and the numerical 
data of Bower [140], and the data of Bradbury [152] are presented in Fig. 7.40. The 
comparison shows good agreement between the computed results and other numerical, 
and experimental results. Comparison between the computed jet centerline velocity decay, 
and other computational results [137, 153, 154], and with experimental data [153] are 
presented in Fig. 7.41. The computed results compare well with experimental data. Data 
presented in Fig. 7.41, were obtained from Jalamani [137]. The discrepancy between 
numerical data are mainly because of difference in inlet jet conditions. 
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V Velocity 


Figure 7.38 Wall Jet Tune- Average- V- Velocity Profiles, at Various Radial Positions, 
Symbol 0 Represents the Spacing of the Grid in the Normal Direction to the Ground Plane 




Static Pressure 


Figure 7.39 Wall Jet Time-Averaged-Static-Pressure Profiles, at Various Radial Positions 
Symbol 0 Represents the Spacing of the Grid in the Normal Direction to the Ground Plane 
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Figure 7.40 Comparison of Ground-Plane Pressure 
Variation for Jet Impinging on a Ground Plane 
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Center Line Velocity Decay, V/V. 
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Figure 7.41 Comparison of Jet Centerline Velocity Decay 
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Jet Impinging on a Ground Plane in Presence of Crossflow When a jet impinges on 
a ground plane a wall jet is formed. The wall jet is symmetric and flows radially outward 
from the zone of impingement, as discussed in the previous section. In the presence of a 
crossflow, the wall jet flowing radially outward is opposed by the crossflow (free stream) 
and rolls up into a horseshoe-ground vortex as shown in Fig. 7 .4. The flow field of the 
impinging jet in a cross flow has been investigated extensively in recent years, because 
of its relevance to the flow field around V/STOL aircraft. 

In the present study, the flow of a jet impinging on a ground plane in the presence of a 
cross flow was simulated numerically. A partial view of the grid used in the computations 
is shown in Fig. 7.42. The grid dimensions were 66x66x33 in the normal, radial, and 
circumferential directions, respectively. The jet Mach number was 0.5 and the Reynolds 
number was 100,000, based on the jet diameter. The jet exit was 3 jet diameters away 
from the ground plane, and the plane extended radially outwards 7.5 jet diameters. The 
ratio of the crossflow velocity to the jet velocity was 0.1 = 0.1). The flow properties 

were normalized with respect to the jet inlet conditions. 

The numerical simulation started by allowing the crossflow to propagate in the 
computational domain until a steady-state solution was obtained. Thus, allowing a 
boundary layer to form on the ground plane. To enhance the rate of convergence of 
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Figure 7.42 Partial View of the Grid Utilized to Compute 
Jet Impinging on a Ground Plane in Presence of Crossflow. 
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the cross flow to steady-state, local time-stepping, residual smoothing and multigrid 
were employed. The next step in the simulation was to release the jet flow into the 
computational domain. Once the jet was initiated into the computational domain, it was 
allowed to interact immediately with the crossflow, and the accelerating techniques; local 
time-stepping, residual smoothing and multignd were turned off for this unsteady flow. 
The solution was advanced in time using a non-dimensional global time step. The time 
step is based on the stability criteria of the scheme. Snap shots, at different time steps, of 
the velocity vectors and the corresponding streamlines are shown in Figs. 7.43 and 7.44. 

The jet propagates into the computational domain and interacts immediately with 
the crossflow, as shown in Figs. 7.43a and 7.44a. The crossflow momentum attempts to 
divert the incoming jet, but due to a lack of sufficient momentum, it fails = 0.1) 
to cause a significant diversion. At the zone of interaction between the free jet and 
the crossflow, a primary vortex starts to form. The size of the vortex increases as it 
is convected toward the ground plane with the jet, as shown in Figs. 7.43b and 7.44b. 
The free jet impinges on the ground plane and is deflected as a wall jet as shown in 
Figs. 7.43c and 7.44c. Subsequently, the stagnation region creates a favorable pressure 
gradient which causes the wall jet to depart from the stagnation region. The wall jet, 
flowing radially outward, is opposed on one side by the crossflow (free stream), and rolls 
up into a horseshoe-ground vortex as shown in Figs. 7.43d and 7.44d. The size of the 
ground vortex increases as it moves radially outward. Figs. 7.43e-h and 7.44e-h. The 
ground vortex forms at a radial location where a momentum balance is reached between 
the wall and crossflow. Fluid mass, from the wall jet (impinging zone) and crossflow, 
accumulates inside the vortex until it cannot sustain itself, and a vortex breakdown occurs. 
A new vortex forms and the process repeats itself. The visualization study conducted by 
Cimbala et al. [147] supports the above description of the flow field of a jet impinging 
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on a ground plane in presence of cross flow. 

A comparison between the computed pressure coefficient distributions, measured 
along the jet centerline, and the numerical results (laminar and turbulent) of van Dalsem 
[142] and the experimental results of Stewart [141] is shown in Fig. 7.45. The present 
results compare fairly well with the experimental data [141], The developed algorithm 
computed the location of the ground vortex accurately. Both the laminar and turbulent 
results of van Dalsem predicted that the ground-vortex was located further upstream than 
the present numerical results, and experimental results. The core of the ground vortex 
corresponds to the point with the minimum C p value. The present numerical simulation 
captured the large scale unsteady phenomena of a jet impinging on a ground plane with 
crossflow, and predicted the location of the ground vortex accurately . Thus the predictive 
capability of the developed algorithm has been demonstrated by accurately resolving the 
complicated fluid dynamics phenomena associated with such a basic, but complex flow 
field. 
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Figure 7.43 Velocity Vectors of the Developing Jet Flow Field in 
Mjet =0.5, Moo =0.1, Pjetl Poo = 1 » Pjet/ Poo = 
Tfi- - l, ^ = 3, Re = lxio 5 (Continued . . . 
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Figure 7.43 Velocity Vectors of the Developing Jet Flow Field in a Crossflow; 
Mjet = 0.5, Moo = 0.1, Pjet/Poo = 1. 

Pj't/Poo = 1. fe = 1. 5 = 3 * Re = 1X105 
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(b) 



(c) 



Moo = 0.1, Pjet/Poc = 1. fe 1 = 1. i = 3, Re = lxio 5 (Continued . . . ) 
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Figure 7.44 Streamlines of the Developing Jet Flow Field in a Crossflow 
M,„= 0.5, =0.1. P,e,/Poo = 1. fe =1-77 =3. fe = lxl0 3 
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Figure 7.45 Comparison of Centerline C p Distributions for a Jet 
Impinging on a Ground Plane with Crossflow; M JC « = 0.5, 

Moo =0.1, P jt tlP x = Ulg = U%=3. He* IxlO 5 . 
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CHAPTER 8 
CONCLUSIONS 


The main objective of the present work was to develop a high-resolution-explicit- 
multi-block algorithm, suitable for efficient computation of three-dimensional, time- 
dependent Euler and Navier-Stokes equations. The developed algorithm employed a finite 
volume approach, and used MUSCL-type differencing to obtain the state variables at the 
cell interface. The variable interpolations were written in the «> scheme formulation. 
Two state-of-the-art techniques were available to construct the inviscid fluxes: Roe s 
flux-difference splitting, and van Leer’s flux-vector splitting. The viscous terms were 
discretized using a second-order central-difference operator. 

The present study investigated two classes of explicit time integration (two-stage 
predictor -corrector schemes, and multistage time-stepping schemes) for solving the com- 
pressible inviscid/viscous flows. The 1-2 Predictor-Corrector Scheme was very effective 
for predicting inviscid flows, but did not perform well with multigrid acceleration tech- 
niques, due to insufficient damping of the high frequency components of the error, over a 
wide range of CFL numbers. The 2-2 Predictor-Corrector Scheme provided much better 
damping for the high frequency components of the error, but the maximum allowable 
CFL number was only one. However, the 2-2 Predictor-Corrector Scheme combined well 
with multigrid acceleration technique and the implicit residual smoothing. 

The standard coefficients of the modified Runge-Kutta method have been modified 
successfully to achieve better performance with upwind differencing. A technique has 
been developed to optimize the coefficients for good high frequency damping at relatively 
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high CFL numbers. The coefficients were optimized using the stability and damping 
factor analysis of the linear wave equation. The optimization was carried out for two- 
, three-, and four-stage schemes. For each scheme the coefficients were optimized for 
four spatial upwind operators: first-order, second-order upwind (k = -1), Fromm second- 
order upwind biased ( k = 0), and third-order upwind biased (* = 1/3), operators. The 
coefficients for a total of twelve schemes were optimized. The multistage schemes 
with the optimized set of coefficients were tested on a number of three-dimensional 
in viscid/viscous flows. In general, the optimum CFL number agreed surprisingly well 
with the stability analysis study, especially for the in viscid cases. 

The explicit upwind schemes, especially for viscous flows, were not effective in 
solving the compressible flow equations, when no additional convergence acceleration 
techniques were incorporated. The main draw back was the restriction on the time step. 
Coupling the explicit schemes with convergence accelerating techniques is essential. The 
accelerating techniques: local time-stepping, implicit residual smoothing, and multigrid 
procedure, proved to be effective in accelerating the rate of convergence to steady-state. 

The implicit residual smoothing extended the stability range of the explicit scheme, 
which simultaneously allowed the use of higher CFL numbers. The smoothing operator 
was applied after every stage of the multistage stepping scheme. Not only did the 
smoothing operator extend the stability range of the explicit scheme, but also provided 
better damping of the high frequency component of the error, which is necessary for 
multigrid acceleration techniques to be effective. The adaptive (variable coefficient) 
implicit residual smoother proved to be successful in accelerating the rate of convergence 
on highly stretched grids with high aspect ratios. 

When using multigrid, the improved damping properties are more important than 
a slight increase in the maximum allowable CFL number. Good damping of the high 


162 



frequency component of the error (close to 0 = x), usually occurred at CFL numbers 
lower than the maximum allowable. Multigrid acceleration techniques enhanced the 
rate of convergence, but true multigrid convergence was not achieved (to the best of 
the author’s knowledge, because true multigrid performance has not been achieved in 
any compressible code). Multigrid performed well with the multistage time-stepping 
techniques and the 2-2 Predictor-Corrector Scheme. It did not perform well with the 
1-2 Predictor-Corrector Scheme, due to insufficient damping of the high frequency 
component of error over a wide range of CFL numbers. 

The developed algorithm was implemented successfully in a multi-block code. The 
multi-block structure provides complete topological and geometric flexibility. The only 
requirement is C° continuity of the grid across the block interface. The solution domain 
was divided into multiple blocks, and the grid for each block was then generated. 
The type of boundary conditions used on each of the blocks was provided through 
an input file at run time, rather than the traditional method of hard-coding the different 
boun dar y conditions in the source code. Non-homogenous boundary conditions per block 
face allowed for more geometric flexibility, and simplified handling the boundaries for 
complex three-dimensional configurations. The application of the developed explicit 
upwind schemes to realistic three-dimensional configurations of significant geometric 
complexity is virtually impossible without the use of multi-block capability. 


The developed algorithm was validated on a number of diverse three-dimensional 
test cases with increasingly complex flow characteristics: (1) supersonic corner flow, (2) 
supersonic plume flow, (3) laminar and turbulent flow over a flat plate, (4) ONERA M6 
wing, and (5) the unsteady flow of a jet impinging on a ground plane (with and without 
cross flow). The test cases were selected to validate certain aspects of the developed 
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algorithm. All of the test cases presented in this study have been computed with the 
same source code; only inputs to the program have been changed. 

Supersonic inviscid corner flow was calculated to verify the multi-block structure. 
It was also used as the bench mark test case for all the explicit time integration 
schemes. Detailed examination of a quasi-two-dimesional jet exhaust plume demonstrated 
that van Leer’s flux-vector splitting, and Roe’s flux-difference splitting gave solutions 
of comparable quality on grids not aligned with shocks and contact discontinuities. 
Extrapolating the conservative variables across slip lines resulted in large over/under 
shoots for Roe’s scheme. The primitive variable extrapolation was found to render a 
smoother solution. 

Laminar and turbulent flow over a flat plate was computed to verify the correct 
implementation of the viscous terms and the Baldwin-Lomax algebraic turbulent model 
in the developed algorithm. The computations demonstrated the effectiveness of the 
adaptive implicit residual smoothing and multigrid procedure in accelerating convergence 
to steady- state. 

The ONER A M6 wing was computed for three different cases: a subcritical test 
case (Moo = 0.699, a = 3.06°), a supercritical case with attached flow (Moo = 0.84, 
q = 3.06°), and a supercritical case with separated flow (Moo = 0.84, a = 6.06°). For all 
three cases, the Reynolds number was 1 1.7 x 10 6 /unit based on the free stream conditions 
and the mean aerodynamic dimension. Good agrement between the computed results and 
other numerical and experimental results was achieved. The developed algorithm cannot 
resolve separated turbulent flows accurately, because of the inadequacy of the Baldwin- 
Lomax algebraic-turbulence model to resolve separated flows. A better turbulence model 
needs to be incorporated if turbulent separated flows are to be investigated accurately. It 
should be noted that turbulence modeling was not the emphasis in this study. 
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The unsteady viscous flow resulting from a jet impinging on a ground plane, with and 
without cross flow, was successfully computed. The different flow fields that identify 
this type of flow, namely the free jet, the stagnation region, the wall jet, the ground 
vortex, and the moving separation zone, were predicted accurately. The jet impingement 
problem demonstrated the capability of the developed scheme to perform time-dependent 
calculations. For unsteady flow computations, if the maximum allowable time step is 
of the same order as the time scale of the physical problem, then explicit schemes will 
best suit the calculations. 

A state-of-the-art computational tool has been developed that is capable of comput- 
ing the flow field around geometrically complex three-dimensional configurations. The 
method has been tested on a number of diverse three-dimensional test cases with increas- 
ingly complex flow characteristics. The developed algorithm is capable of accurately 
resolving steady and unsteady compressible flow fields. To further enhance the perfor- 
mance of the developed algorithm the following recommendations are suggested: 

The damping of the explicit time-stepping schemes on highly stretched grids with 
cells that have extremely high aspect ratios should be investigated. A stability analysis 
study should be conducted on the two-dimensional, scalar wave equation and the diffusion 
equation to investigate the effect of the gnd aspect ratio on the damping characteristics 
of the multistage schemes. This study should take into account the effect of implicit 
residual smoothing and multigrid procedure. 

Multigrid acceleration techniques performed well for inviscid non-stretched grids, 
but the convergence rate deteriorated on viscous grids. One possible remedy lies in 
using the technique of semi-coarsening in the direction normal to the wall, to improve 
the non-uniformity on the coarse-grids. Other approaches have to be combined with the 
multigrid method to damp the high frequency-component of error at all grid levels. An 
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attempt to develop a mathematical operator to annihilate the high frequency component 
of the error should be investigated. 

A joint analytical and numerical study should be conducted to develop higher order 
accurate, non-oscillatory schemes. A stability study should be performed to investigate 
the damping characteristics of the higher order schemes. 

Explicit schemes are amenable for implementation on massively parallel supercom- 
puters. Work is already in progress to implement the developed algorithm on massively 
data-parallel architecture (Connection Machine CM-2). 

A higher order turbulence model needs to be employed in the developed algorithm, 
to be capable of resolving turbulent separated flows. 
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Appendix A 

Full Navier-Stokes Equations 
in Body-Fitted Coordinates 


The time dependent compressible three-dimensional Navier-Stokes equations in gen- 
eral curvilinear coordinates, written in strong conservation form (neglecting the body 
forces and the external heat sources are [6] 

d_Q d{F-F t ,} d{G-G v ) d{H-H v ) Q (A1) 

dt + di d( 


where Q is the vector of dependent variables and is given by ; 
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The inviscid fluxes F, G, and H are function from the state vector Q . They are given by : 
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where U, V , and W are the contravarient velocity components in the G rj and (. They 


are given by , 
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The viscous flux, F v is given by F v = {F Vl , F V2 , F V3 , F Vt ,F Vs }, where. 
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The viscous flux, G v is given by G v = { G , G Vi , G „ 3 ,G Vt , G Vb } where. 
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u di^Vx + ‘yVy + Zzilz) + u v {%vl + Vy + vl) + 


, M re f ft 

* v 2 — 

i\ T e f 


UclfCz 7 ?! + CyVy + CzVz) + 
VtitiVy ~ jty'h) + v y^ T h^ly + ^(Cx»?y - %(yVx) + 
W({ZzT]z ~ f ZzVx) + w t l^VxVz + ^(Cl^z — |C zVx) . 


(A.9) 


(A. 10) 


(A.ll) 
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G 


1'3 


G v 4 


M re f ft 
Rref 


M re f ft 

Ft r e f 


udtyVx ~ \ixVy) + U V3W X + u c(C y^x ~ 3^9) + 

v ( i {it T ix + f£t rfy + ZzVz) + v y{ r lx + 3^y + ^z)^" 

V((C xT]x + fCy^y + CzVz) + 
Wt(tyT) z - J izVy) + Wv^VyVz + «-’C {CyVz ~ f CzVy) - 

'utitzVx - § ZxVz) + u AVzVx + u c {(zVx ~ | CzVz) + 
v t(ZzVy ~ § iyVz) + v v \i izVy + v c(^ zT >y ~ 3^*) + 

Wz(£, x T) x + i y T]y + \izVz) + w y {jix + % + 3^) + 

W(((iT) x + (yTjy + j(zVz) . 


(A. 12) 


(A. 13) 


Cl Vi — 


MrefP 


R 


ref 


u z 


li^x^lx + ^yly iz^z) + v i^y 7 1 x 3^1 Vy) J r 

u>{£zVx ~ f ixVz) J 


+ 


l„[u{hx + O + 3 V Wx + \wizVz] + 

u c 


u(| CxVx + ^yVy + (zVz) + v {Cy r lx 3(3:^) + 

w{CzV x - ICxVz) J 


u(^j-'/v — f^y 7 /i) + v {tx 7 1x + 3 ZyVy + (,zVz) + 
"<[ w{ZzVy~UyVz)\ 

>„ [itll/,1?, + I'ijVl + 02 ) + l wr l9Vz] + 

u(G'/ V ~ fC yVx) + v ((xVx + jCyVy + C*»7*) + 

w{(zVy ~ iCyVz) 

u{£,x 7 h — fs zVx) + v ((y 1 lz — §£::%) + + 

w (£i 7 1i + £yVy 3 Zzlz) _ 

XV,,[\uj) : 1] x + \vilzVy + w ih 2 + ^)] + 


+ 

+ 


+ 


w t 


W C 


u{OVz - fC zVx) + v( Cy*iz - | CzVy) + 
w((xVx + C yly + fC zVz) 


+ 


_cr[Tz(Z z T h + ZyVy + ZzVz) + T„e 2 + T c ( ( x Vx + Cy% + 


(A. 14) 
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The viscous flux H v is given by, H v — {//«,, Hvi i , H Vs } where. 


H Vl = [ 0 ] 


(A. 15) 


_ Mrc fH 

R-re f 


u {(fGG + GC y + GG) + U v (jT) x ( x + T}y( y + T)z(z) + 

u c(|Ci + Cy + Cz) + 

V*U*G - §GG) + U|»(»/xC» - §»?wG) + «ClGG + 
u^(GG — §GG) + wt){vxCz — fVzCx) + ™cjGG 


A/ r e//^ 

/£ rf y* 


u*(GG - §GG) + «n(»?»G - i^G) + u ciGG+ 
t'{((xG + §GG + GG) + u ?j(^iG + f^yG + r ?xG)+ 

u c(Ci + fCj + G) + 


u>*(GG - |GG) + w viv y Cz - | izCy) + u 'cfGG 


(A. 17) 




Mrt f H 

Rref 


u ( (GG _ fGG) + u,,(t) z ( x — 3 ^xG) + MC 3GG+ 
v*(GG - fGG) + t’„(»jxC* _ I 7 ?/*) + v cfGG+ 
U)(((iG + GG "b fGG) + w r}{ T lx G "b tfyCy + 3 7^G ) “b 


(A. 18) 


™<(C* + Cj + fG) J 
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H Vi = 


Mre f M 
R-ref 


H 


+ 

+ 


(f^xCz + Zy(y + £zCz) + v ((y(z ~ ;?£zCy) + 

w (£z( x — f^iCz) J 

' u(±T) x ( x .+ T)y(y + Vz(z) + v (Vy(x - f VzCy) + 

u "[ w(r, s Cz-b^)\ 

U C [ u (iC + , 'T >2 ) + I^CyCr + J^CzCx] + 

u{Zx(y — 3 Zytx) + u(6r(, z + f^yCy + £zCz) + 

w {tz(y — 3 £yCz) J 

u i 7 UCy ~ f%Cr) + v { T h(x + f 9yCy + »?*Cz) + + 

w{i]z^y ~ s^yCx) . 

u c[3 u C»Cr + u (Ky +V’ 2 ) + 3 u, CzCy] + 


v t 


+ 


W, 


l^xG — f^zCx) + v {vy(z ~~ f*?*Cx) + 
w {jlx Cx + lyCy "t" ^zCz) J 


+ 


U ’C [3 “GCr + jK-'Cy + MK* + ( r ,2 )] + 

.vfcitxCx + £yCy + 6Cz) + T v (Vx(x + riy(y + VzCz) + T^ 2 ] 


(A. 19) 
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Appendix B 

Derivation of the Reynolds-Averaged 
Navier-Stokes Equations 


The Reynolds-Averaged Navier-Stokes equations are derived from the Navier-Stokes 
equations by decomposing the randomly changing flow variables into a mean and a 
fluctuating components. For compressible flow triple correlations involving density 
appears in the equation. Favre [98] suggested a mass weighted decomposition for 
compressible flows to avoid the triple correlation involving density. The following 
formulations were used to decompose the flow variables in the Navier-Stokes equations 


2.1. 


7 p£ ~ = f = — H - — 

p ' u p ’ p ’ p 

where u = u + u, T = T + T, H = H + H 
note p = p + p, p — ~p + p 


(B.l) 


where (p + p)f = 0 
but f / 0 

By substituting these quantities into the governing equations and averaging in time , we 
get the continuity equation 


^ + ^ = 0.0 

dt dx t 


(B.2) 


the momentum equations 

d[pui\ d{pu t u } } _ _ dp_ + - pujUj } + _ df 

dt dxj 3x x 3xj dx t 


(B.3) 
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the energy equation 


d^pu t H + pu\Hi — ujTij — TijU t 
dt + dx, 


it can be shown that 


0.0 


(B.4) 


H = h + I\ + K 


(B.5) 


and 


H = h + u,u, + A — h 


(B.6) 


with 


-pK = = pk, and k= (B.7) 

where K is the average of the turbulent Kinetic Energy of the Turbulent fluctuations 
which by definition is K. Hence 


u,u. 


pu 


iH = piixh - r£ n Um + puil< 


(B.8) 


where 


T,m = -PH, Urn 


(B.9) 


and is defined as the Reynolds stress tensor. Substituting into the energy equation gives 


»{e} 


+ ^-\pu,H + pu,h - k'T~- - T,jUj - r,* u m - T l} Uj + pu,I\ 1- = 0.0 (B.10) 
dt dx, \ dx, ) 

A host of terms involving the fluctuations evolved and need to be modeled in order to 

solve the governing equations. These terms are 


dT 


pu,h. T t jUj, T,jU , , pu, K 


(B.l 1) 


pu,h is generally considered the turbulent heat flux, and is modeled as 


pu,h = —kjVT 


(B.12) 
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where Jb T is the coefficient of turbulent thermal conductivity. r t j will be modeled using 
Boussinesq’s approximation [98]: 

-IpKSij (B.13) 

is the turbulent eddy viscosity and is related to k x by 


T-j = -puxUj = pi 


’ 3uj du, 2 du m 

dx{ dxj 3 dx. 


o 


= 


CpP T 

P rT 


(B.14) 


p T T is the turbulent Prandtl number and is equal to 0.9 . r t j and r l} are combined 
together to give 


= (p + pt) 


duj du t 2 du r 


+ ~ 


[d.r, dxj 'idx r 




- 1~PK8 XJ 


(B.15) 


In the present work t X jU } and pu x I\ are assumed to be small and are neglected. They 
are calculated in a higher order scheme. 

Applying the same procedure to the equation of state we get 


P 


h - i) 



pu t u, 

—- pk 


(B.16) 


Thus, the Reynolds-Averaged Navier-Stokes equations in dimensional Cartesian form 


become; 


where 


dq_ d(h-h v ) 

dt c)x dy dz 


< 7 


p u 

< ~pv > 

pw 

. E . 


(B. 17) 


(B.18) 
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p u 

— ~2 t — 

pu- + p 


/=< j put' >, /„ = (/x + /iT)^ 

p uu? 

, (c + p)u > 


uf 


odu _ dr 
“ dx dy 


o.o_ 

dw 

dz 


+ 


2 pK 

3(#i+MT> 


O du 

2 3F “ 


dr 

3F 


dr i du 
dx^ ' dy 
dw i du 

~ 5 x \Tfz 


dw 

1 z 


+ 


~( dv 

V 3 " 


+ £) 


+ 


w 


(fr + li) + u 3 (m+/i T ) + _ 

(*+^Xl dT 

(p+jxt) # 


(b.19) 


= < 


pv 

p uv 

— — -2 , — 
p l’ + p 

Jvw 

{ (e + pit 


>, 9 v = (p + Mt)s 


J>.0 _ 

du i dv 

+ Tfx 

l\‘~)dv du dw 

’ '^Hy Hx ~$z 


+ 


2pK 


3 (m+/jt) 


dv I 
dz i " 


du , d^\ , ~2 o9t> _ du _ dt£ 

U ^+dFj +U 3[ Z % dx dz 

2 pK 


dw 

d t 

dv 

dy 


w 


( dv i duA 

^dz 


+ U 


+ 


(k+k T ) dT 


(m+mt) d 


!b.20) 





h v = in + Pt)\ 


^0.0 _ 
du i dw 
d £ du 
dv 1 dw 
dz ' dy 

2 odir _ du _ dr i 'ZpK 

3 “ dz dz dy ' 3(/i + /iT) > 



dr . dir A , 
dz + dy 


dv 

dy 


i ~ 2 p K | 

U 3 (fi+flT ) 


(tthidT 
lM+ “ Tl %.21) 


the superscript .” has been removed for clarity, and the lower case letters 
q. /. f v , 9 - 9 v h. and h v are used to identify these vectors are functions from 


the averaged variables. 
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In this study, the Baldwin-Lomax, algebraic-turbulence model was selected which 
model the eddy viscosity. Baldwin-Lomax model implements a simple algebraic expres- 
sion for the eddy viscosity. The model is easy to implement and provide reasonable 
results for attached flows. For separated flows or for the purpose of investigating the 
turbulent behavior of fluid motion, it is recommended that a higher order models, such 
as the two equation models I\ - e, K - w, and K - r [155] or solve for the Reynolds 
stresses [6]. The Baldwin-Lomax model is implemented with the understanding of the 
limitations and restrictions of the model. Higher order models should be included in the 
code to give the generality and accuracy in solving turbulent problems. 

There is no mechanism in an algebraic-turbulence model that can account for the 
turbulent kinetic energy,/\ . Thus it was dropped from the governing set of equations as 
well as from the equation of state. The Reynolds-Averaged Navier-Stokes equations in 
non-dimensional Cartesian form takes the following form ; 
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G={ 


pv 

puv 

pv 2 + p 
pvw 
{ {E + p)v ) 


G v — 


M 


ref 


0.0 

du i dv_ 
^ dx 


R 


ref 


(M + MT ){ 


dy ' dx 

2 ( o dx: du _ div \ 

3 dy dx dz ) 


dv I dw 
dz ' dy 


( du , dv \ . 2 / 9 _ du I 

u [dt + 37 J + ^3 V“^y ^ d: ) + 

( dv , d?i> \ I I 

^ ^. 24 ) 


// 


pw 
puw 
pvw 
pw 2 + p 
{(E + p)w ) 


M 


H v= ‘-pL(p + p T ){ 

Kref 


0.0 

du . du? 

lit Is 

dz ' 


2/^9 a™ _ du. _ dv \ 

3 y* dz dx dy ) 

.. ( du j_ dw \ i 7 ( dv_ I duA 

u ld7 + 37) +^37 + a y j 

| 2 / 9 du* du dv \ i dT 1 

{ w s{ 2 ^-d;-di) +(T dz ) 


(B.25) 


where 


(7 = 


i (f_ + JfL\ 

3(p + pt)\P t PtT ) 


(B.26) 
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Appendix C 
Limiters 


C.1 Minum-Modulus Limiter 

The min-mod limiter is implemented in the present algorithm by rewriting the 
K-scheme formulation, eq. (19), as follows 

+t[(1-«)V 1 + (1 + k )A 1 ] 

,+ * 9 (C.l) 

Q* x = Qi + 1 - 1 [ (i - «)V i+1 + (l + «)A i+1 ] 

where y, and A r are defined as 


A, = minmod( A,,/iy,) 
y, = minmod(SJi, 3A t ) 
At — Qi+i Qi 


(C.2) 


y. = Qi - Qi-\ 

where. A, and y, are the backward and forward differences, and 3 is the compression 


parameter given by 


3 = 


:5 - k 

1 - K 


(C.3) 


The min-mod operator is 


minmod(x, y) = sign(x) maj-{0.mrn|x sign{y},y «0 n{x}|} (C.4) 
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C.2 Vanalbda Limiter 


The Vanalbda differentiable limiter modifies the k — scheme formulation, eq. (19), 
as follows 

0‘ , =Q , +^[(1-kS 1 )V 1 + (1+kS,)A,] 

1+ = 4_ (C.5) 

o* , = <3,+, - ^r 1 ! (1 - -cSi+OV*. + (1 + K5 j+ ,)A i+ , ] 

where 

2 A, y, + f (C.6) 

vr + + < 

and e is a small number to prevent the division by zero in the region of zero gradients. 
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Appendix D 
Multigrid Cycles 


Three types of cycles have been employed in this study, V-Cycle, W-Cycle, and the 
Full Multigrid Cycle (FMG). Basically, a standard V-cycle can be broken into halves. 
The first half is the restriction part of the cycle going from the fine-grid through the 
coarser grids down to the coarsest grid. The second half is the prolongation part of the 
cycle going from the coarsest grid up to the finest grid. An example is shown in Fig. D.l 
for a four level multigrid. The circles indicate when iterations are performed on the given 



grid level, and the lines between grid levels indicate either a restriction or prolongation 
operation. Notice that the circle for the fine-grid at the beginning of the cycle is omitted 
since the iterations on the fine-grid are performed at the end of the prolongation section. 
This ensures that the last operations in a multigrid cycle include updates on the fine-grid. 
It is often necessary to perform more than one iteration on a given grid to get the required 
smoothness in the error for multigrid to work. 

A W-cycle can be thought of as consisting of several components which are similar 
to V-cycles but with different varying ‘coarsest’ and ‘finest grids This idea is shown 
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in Fig. D.2 where a W-cycle is graphically expanded to show its ‘legs’. This allows a 
simple coding modification to the V-cycle program to allow W-cycles. 



Figure D.2 W-Cycle. 

The Full Multigrid Cycles (FMG) are used to get a good initial approximation on 
the fine-grid, which can be used as the starting solution for a V- or W-cycle. Figure D.3 
shows the schematic of a FMG V-cycle, while the schematic of a FMG W-cycle is 
shown in Fig. D.4. The basic idea of a FMG cycle is to iterate first on the coarsest 
levels, assuming the finer level to be the solution level. The solution is then prolongated 
to the next finer level and the problem is solved again on these three levels. The process 
is repeated until the finest level is reached. The FMG cycles are inexpensive since we 
are solving the problem on a coarse-grid. By the time we include the finest grid in the 
FMG procedure and start to apply the regular V- or W-cycle, the global features of the 
solution has already been developed. 
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